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SECTION  1 


PROJECT  OVERVIEW 

1.1  SUMMARY  OP  OBJECTIVES  AND  METHODOLOGY 

The  response  of  saturated  soils  and  rocks  to  high  intensity  blast 
loadings  is  a  subject  of  intense  concern  to  the  United  States  Air  Force. 
Previous  studies  (Blouin  and  Kirn,  1983;  Blouin  and  Kim,  1984a;  Blouin  and 
Shinn,  1983)  clearly  defined  the  role  of  blast  induced  liquefaction  on  the 
cratering  processes  from  both  nuclear  and  high  explosive  charges  and  the 
subsequent  flow  and  consolidation  of  liquefied  material  surrounding  and 
beneath  the  crater.  These  studies  have  been  used  by  the  Air  Force  and  the 
Defense  Nuclear  Agency  (ONA)  to  help  explain  the  anomalously  large,  shallow 
craters  which  were  observed  in  the  saturated  geologies  of  the  Pacific  Proving 
Grounds  where  all  U.S.  high  yield  nuclear  surface  tests  were  fielded. 

In  addition,  prediction  of  the  ground  shock  transmitted  by  a  surface  or 
shallow  buried  explosion  to  nearby  surface  or  deeply  buried  hardened  struc¬ 
tures  is  a  key  element  in  the  design  and  evaluation  of  such  structures.  In 
order  to  make  reliable  ground  shock  predictions  in  saturated  geologic  media, 
use  of  multiphase  material  models  and  numerical  codes  is  needed. 

This  report  summarizes  results  of  a  three  year  combined  experimental  and 
theoretical  study  of  the  fundamental  behavior  of  multiphase  porous  materials 
subjected  to  high  intensity  dynamic  loadings.  The  overall  and  specific  objec¬ 
tives  of  this  study  were; 

1.  To  design  and  conduct  laboratory  experiments  which  identify  key 
response  mechanisms  and  measure  the  behavior  of  saturated  porous 
materials.  Specific  areas  of  concentration  under  this  task  included: 

a)  design,  fabrication  and  utilization  of  a  device  to  measure  fluid 
friction  and  inertial  flow  resistance  through  ducts  and  soil  and 
rock  samples  at  high  pore  pressure  gradients  including  flow  in  the 
laminar,  transitional  and  turbulent  flow  regimes; 
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b)  measurements  of  the  compressibility  of  soil  and  rock  grains, 
including  the  influence  of  occluded  porosity  on  the  grain 
compressibility; 

c)  documentation  of  the  liquefaction  process  in  both  soils  and  porous 
rocks  under  uniaxial  strain  loadings  similar  to  those  produced  by 
explosive  loadings;  measurement  of  the  amount  of  consolidation  of 
the  soil  and  rock  specimens  following  liquefaction; 

d)  experimentally  determine  the  amount  and  severity  of  grain  breakage 
as  a  function  of  both  strain/stress  path  and  stress  magnitudes; 
relate  the  microscopic  grain  crushing  response  to  the  degree  of 
shear  and  compression  generated  during  the  various  types  of 
loadings;  and 

e)  experimentally  determine  the  drained  skeleton  response  and  the 
undrained  total  stress,  pore  pressure  and  effective  stress  response 
under  various  simple  and  complex  loading  conditions. 

2,  The  development  of  advanced  fully  coupled  material  models  which 
accurately  represent  the  dynamic  multiphase  response  behavior  of 
saturated  and  partially  saturated  soils  and  rocks.  Specific  areas  of 
concentration  under  this  task  included; 

a)  formulation  of  a  pore  fluid  flow  relationship  which  describes  the 
pore  fluid  flow  resistance  under  dynamic  loading  conditions  in  both 
the  laminar  and  turbulent  regimes;  and 

b)  development  of  a  procedure  to  model  the  influence  of  occluded 
porosity  on  the  response  of  saturated  and  partially  Saturated 
soils  and  rocks; 

3.  formulate  and  implement  both  the  theoretical  models  and  experimental 
results  into  two-phase  and  multiphase  codes.  Specific  accomplishments 
under  this  task  include; 
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a)  development  of  the  numerical  code  NKOCP  which  is  used  to  numeri¬ 
cally  model  the  hydrostatic  and  uniaxial  undrained  response  of 
saturated  porous  materials; 

b)  development  of  the  code  TWAVE  which  provides  rapid  closed  form 
solutions  of  wave  propagation  and  damping  in  saturated  porous 
materials  having  linear  elastic  skeletons,  including  definition  of 
compress ional  waves  of  the  first  and  second  kind; 

c)  development  of  the  numerical  finite  element  code  MPOAP  which  utili¬ 
zes  a  single  point  stress  calculational  technique  and  an  innovative 
method  of  computing  pore  fluid  flow  to  maximize  computational  effi¬ 
ciency;  features  of  MPOAP  include: 

•  nonlinear  fluid  friction  model  for  both  laminar  and  turbulent 
flow; 

•  a  fully  coupled  compressibility  model; 

•  a  number  of  drained  skeleton  material  models  ranging  from  simple 
linear  elastic  to  advanced  elasto-plastic  models  with  strain 
hardening  and  softening; 

•  single,  two  and  three-phase  capability; 

•  static,  dynamic  and  quasi-static  capability  with  pore  fluid  flow; 

.  orai  - 

'  ‘  N. 

d)  development  of  modifications  to  our  existing,  but  less  sophisti¬ 
cated  two-phase  finite  element  code  TPOAPII.  These  modifications 
were  developed  as  a  result  of  the  MPOAP  theoretical  formulations 
and  were  used  in  many  of  the  numerical  studies  under  this  project. 

4.  Utilize  numerical  and  theoretical  codes  to  help  analyze  experimental 
results,  to  predict  laboratory  and  field  behavior,  and  to  conduct 
numerically  based  experiments  and  parameter  studies  which  will  further 
advance  our  understanding  of  complex  multiphase  phenomena.  Specific 
accomplishments  under  this  task  include: 

a)  modeling  of  undrained  laboratory  response  of  saturated  soils  and 
rocks  to  a  variety  of  load-unload  conditions,  including  the 
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modeling  of  liquefaction  in  soils  and  rocks  due  to  single 
hydrostatic  and  uniaxial  strain  loading  cycles; 

b)  numerically  model  wave  propagation  in  saturated  soils  and  rocks  due 
to  dynamic  loadings,  including  definition  of  the  liquefaction  front 
and  motion  of  the  pore  fluid  as  a  function  of  position  behind  the 
wave  front; 

c)  theoretical  calculations  of  wave  propagation,  velocity  and  energy 
damping  as  a  function  of  the  skeleton  permeability  and  excitation 
frequency;  define  the  existence  of  compressional  waves  of  the 
first  and  second  kinds  and  determine  propagation  velocities  and 
damping  differences  between  the  two  wave  types;  and 

d)  numerically  model  the  theoretical  wave  propagation  results  in  *'c" 
above  at  specific  values  of  excitation  frequency  and  permeability 
to  determine  the  role  of  pore  fluid  motion  on  damping  and  wavespeed 
and  to  isolate  the  physical  phenomena  which  govern  propagation  of 
waves  of  the  second  kind. 

An  integrated  experimental/theoretical/calculational  approach  has  been 
used  to  satisfy  the  above  objectives.  In  this  approach,  each  of  these  three 
aspects  of  the  program  is  mutually  supported  by  the  other  two.  For  example, 
in  the  shock  consolidation  laboratory  tests  reported  by  Kim,  Blouin  and 
Timien,  1981,  liquefaction  of  porous  limestone  and  apparent  negative  hystere¬ 
sis  in  the  pore  pressure-volume  strain  curve  were  both  observed.  Numerical 
simulations  of  those  experiments  utilizing  the  NKOCP  code,  developed  from  our 
fully  coupled  two-phase  model,  allowed  us  to  duplicate  the  experimental 
results  and  identify  the  governing  two-phase  phenomena  which  produced  these 
surprising  results.  On  other  occasions,  we  have  identified  theoretical  uncer¬ 
tainties,  such  as  dynamic  pore  fluid  friction,  and  designed  and  analyzed 
experiments  to  define  the  controlling  equations  and  material  property  parame¬ 
ters  used  in  the  theoretical  formulation  and  numerical  implementation.  When 
used  in  these  ways,  the  combined  experimental/theoretical/numerical  approach 
has  been  a  powerful  technique  for  solving  problems  which  might  prove  intrac¬ 
table  to  solution  by  any  single  element  of  these  three. 
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This  report  is  the  third  and  final  report  on  this  study.  There  have  been 
two  annual  reports  by  Kim,  Blouin  and  Timian,  1986  and  1987,  summarizing  deve¬ 
lopments  at  the  end  of  the  first  and  second  years,  respectively.  In  some 
cases,  the  material  presented  in  those  reports  is  complete  and  will  not  be 
repeated  in  this  report.  Rather,  a  short  summary  of  the  significant  findings 
from  these  two  reoorts  is  included  in  Subsection  1.2.  In  other  cases, 
material  presented  in  these  reports  has  been  supplemented  or  finalized  during 
the  last  year.  In  these  cases,  the  final  development  of  these  materials  is 
included  as  a  section  of  this  report,  with  the  previous  work  incorporated  into 
this  final  presentation. 

Key  individuals  responsible  for  this  project  at  ARA  include  the 
following: 

*  Or.  Kwang  Jin  Kim:  Co-Principal  Investigator,  Senior  Engineer  -  theore¬ 
tical  analysis,  numerical  implementation  and  analysis,  experimen¬ 
tal  analysis; 

•  Hr.  Scott  E.  Blouin:  Co-Principal  Investigator,  Principal  Engineer  - 
experiment  design  and  analysis,  theoretical  and  numerical 
analysis; 


•  Or.  Douglas  Kerkle:  Principal  Engineer  -  theoretical  analysis? 

•  Hr.  David  Timian:  Staff  Engineer  -  experiment  design  and  analysis; 

•  Hr.  Oaniel  Chitty:  Senior  Engineer  -  experiment  design  and  analysis; 

•  Hs.  Elizabeth  Smith:  Staff  Engineer  -  numerical  analysis: 

«  Hr.  Lawrence  Kerkle:  Technician  -  numerical  development  and  analysis, 

•  Mr.  Kenneth  McIntosh:  Staff  Technician  -  experiment  design,  execution 
and  analysis; 

•  Hr.  Steven  Quenneville:  Technician  -  experiment  design,  execution  and 
analysis; 
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•  Hr.  Douglas  McIntosh:  Technician  -  experiment  design,  execution  and 
analysis. 

A  number  of  formal  presentations  have  been  given  presenting  work 
accomplishment  under  this  study  or  presenting  results  and  analysis  which  were 
strongly  supported  by  this  study.  Among  these  are: 

1)  "Fundamental  Analysis  of  Wave  Propagation  and  Liquefaction  in  Multiphase 
Porous  Media"  by  Scott  E.  Blouin,  Kwang  J.  Kim,  David  A.  Timian;  Presented 
at  AFOSR  Soil  Mechanics  Seminar,  MIT,  14-15  September  1987. 

2)  "Strength  and  Deformation  Properties  of  Salem  Limestone"  by  Scott  Blouin 
and  Oaniel  Chitty;  Presented  at  ONA  Material  Properties  Working  Group, 
Weidlinger  Associates,  Mew  York,  NY,  19  July  1988. 

3)  "Some  Aspects  of  Two-Phase  Modeling  Investigated  Under  Sponsorship  of 
AFOSR"  by  Kwang  Kim  and  Scott  Blouin;  Presented  at  DMA  Mange  to  Effect  EPW 
Meeting,  SAIC,  Albuquerque.  NH,  17-18  May  1988. 

4)  "Comparisons  and  Calculations  of  Laboratory  Properties  of  Various 
Limestones"  by  Scott  81ouin,  Kwang  Kim  and  James  Drake;  Presented  at  ONA 
Mange  to  Effect  Meeting,  Los  Alamos  National  Laboratory,  20  January  1988, 

5)  "Preliminary  Comparisons  of  Salem  and  ONA  Limestone"  by  Scott  81ouin,  Kwang 
Kim  and  Robert  Walker;  Presented  at  ONA  Range  to  Effect  Meeting,  SAIC, 
Albuquerque,  NM,  17-18  May  1988. 

In  addition  to  the  above  and  other  formal  presentations,  research  results 
from  this  study  have  been  supplied  on  a  continuing  basis  to  universities, 
research  companies,  national  laboratories  and  government  organizations  in  the 
form  of  copies  of  briefings  and  annual  reports,  and  consultations.  Government 
organizations  and  national  laboratories  to  which  we  have  supplied  information 
include: 

Defense  Nuclear  Agency 

Air  Force  Weapons  Laboratory 
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Air  Force  Engineering  and  Services  Laboratory 

Los  Alamos  National  Laboratory 

Sandia  National  Laboratory 

Lawrence  Livermore  National  Laboratory 

1.2  SUMMARY  Of  SIGNIFICANT  FINDINGS  AND  ACCOMPLISHMENTS 

In  Subsections  1,2.1  and  1.2.2,  we  briefly  summarize  the  significant 
findings  and  accomplishments  which  are  fully  described  in  our  previous  annual 
reports  and  which  are  not  updated  and  described  in  Sections  2  through  ?  of 
this  report.  References  to  appropriate  sections  in  the  previous  annual 
reports  are  included.  In  Subsection  1.2.3,  we  briefly  summarize  the  findings 
presented  in  Section  2  through  7  of  this  report. 

1.2.1  Summary  of  Laboratory  Experiments  and  Analysis  from  Previous  Annual 
Reports 

a.  Grain  Compressibility  and  Influence  of  Occluded  Porosity  -  Kim, 

Blouin  and  Timian,  1906  Sections. 

The  compressibility  of  the  solid  grains  is  an  integral  part  of  the 
constitutive  formulations  for  multiphase  response  (see  Section  2).  The 
compressibilities  of  the  mineral  constituents  such  as  quartz  and  calcium  car¬ 
bonate  making  up  most  of  the  soils  of  interest  are  generally  available  in  the 
literature  from  sources  such  as  Bridgman,  1931  and  Simmons  and  Wang,  1971, 
During  a  previous  study  Qlouin  et  a!,,  1904,  observed  that  the  grains  of  the 
carbonate  soils  and  rocks  from  Enewetak  Atoll  contained  a  high  degree  of 
microporosity  within  the  grains  themselves.  The  pervasiveness  of  these  micro¬ 
pores  explains  why  in  situ  densities  of  many  carbonate  soils  and  rocks  are  so 
low.  A  porosity  value  of  SOX  is  typical  of  both  the  uncemented  and  cemented 
sediments  at  Ebewetak.  In  the  case  of  Enewetak  beach  sand,  perhaps  a  third  of 
the  bulk  porosity  is  due  to  the  intragranular  microporosity. 

If  the  micropore  space  in  the  carbonate  sediments  is  not  fully  saturated, 
resultant  undrained  compressibilities  can  be  much  greater  than  those  predicted 
using  the  solid  grain  compressibilities  from  the  literature.  In  essence,  the 
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bulk  moduli  of  the  solid  grains  could  be  greatly  reduced  by  the  presence  of 
the  unsaturated  micropores,  leading  to  a  substantially  more  compressible  bulk 
mixture.  Under  a  high  stress  load-unload  cycle,  such  a  mixture  would  exhibit 
permanent  compaction  and  high  energy  absorption,  neither  of  which  are  charac¬ 
teristic  of  fully  saturated  materials.  Thus,  it  is  important  to  characterize 
the  grain  compressibility  of  materials  with  microporosity  to  insure  that 
appropriate  relationships  for  grain  compressibility  are  available. 

A  laboratory  high  pressure  vessel  was  modified  to  enable  us  to  measure 
the  compressibility  of  known  volumes  of  fluid  (kerosene)  within  the  vessel. 

By  mixing  predetermined  amounts  of  soil  or  rock  grains  in  the  fluid,  the 
compressibility  of  the  grains  can  be  determined  by  measuring  the  compressibi¬ 
lity  of  the  fluid-grain  mixture.  There  is  no  effective  stress  applied  to  the 
grains  in  this  test;  the  only  stress  acting  on  them  is  the  fluid  pressure. 

The  grain  bulk  moduli  of  the  various  soils  and  rocks  examined  in  this 
suite  of  tests  are  summarized  in  Table  The  grains  all  exhibited  essen¬ 

tially  linear  elastic  response  to  the  peak  applied  pressure  of  about  5  kb. 
Quartz  sand,  steel  ball  bearings  and  solid  limestone  produced  bulk  moduli 
which  agree  well  with  published  values.  Most  of  the  soils  and  rock  grains  had 
bulk  moduli  of  about  6.0  X  105  psi,  about  55%  of  the  modulus  of  the  solid 
limestone.  While  this  is  apparently  a  significant  reduction,  in  reality  it 
will  have  little  influence  on  the  response  of  the  Enewetak  materials  because 
the  water  making  up  50%  of  the  sample  is  more  than  an  order  of  magnitude  more 
compressible  than  the  grains.  Thus,  it  was  concluded  that  the  micropores  are 
probably  mostly  saturated  in  these  materials  and  result  in  a  20%  to  45%  reduc¬ 
tion  in  grain  modulus. 

b.  Liquefaction  of  Saturated  Soil  and  Rock  Under  Undrained  Uniaxial 
Strain  Loadings  -  Kim,  Blouin  and  Tiraian,  1906  Section  3,  1987 
Section  7. 

A  series  of  uniaxial  strain  load-unload  consolidation  tests  was  run  on 
porous  limestone  and  soils  from  Enewetak  Atoll.  The  test  specimens  were  sub¬ 
jected  to  an  undrained  load-unload  cycle  in  a  high  pressure  oedometer  during 
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which  pore  pressure  was  monitored.  Following  unloading,  the  pore  pressure  was 
allowed  to  drain  until  the  samples  had  consolidated  to  their  original  effec¬ 
tive  vertical  stress  at  the  start  of  each  test. 

Figure  1.1,  parts  a,  b  and  c,  show  a  typical  data  set  on  a  saturated  low 
strength  limestone  having  a  porosity  of  42%.  Part  a  shows  the  total  stress 
response,  part  b  the  pore  pressure  response  and  part  c  the  effective  stress 
response.  Point  "a"  at  the  start  of  the  undrained  loading  shows  the  initial 
conditions,  i.e.  effective  stress  of  about  300  psi,  pore  pressure  of  about 
300  psi  and  total  stress  of  600  psi.  The  sample  is  loaded  undrained  to  point 
"c".  At  point  "b",  early  in  the  undrained  loading,  the  cementation  in  the 
skeleton  breaks  down.  This  is  clearly  shown  by  the  effective  stress  response 
of  Figure  1.1c.  Prior  to  the  breakdown  of  cementation  more  than  half  of  the 
total  applied  stress  is  carried  by  the  limestone  skeleton;  but  following  the 
cementation  breakdown  only  5%  of  the  additional  applied  stress  is  carried  by 
the  skeleton.  The  balance  of  the  total  applied  stress  is  carried  by  the  pore 
water. 

From  point  "c",  the  undrained  sample  is  unloaded  to  a  state  of  zero 
stress  at  "e".  As  shown  in  Figure  1.1,  the  skeleton  is  strongly  hysteretic 
and  unloads  very  rapidly  until  at  "d"  where  the  effective  stress  drops  to  zero 
and  a  state  of  liquefaction  is  achieved.  Note  that  at  point  "d"  the  overall 
volume  strain  is  still  1,4%  and  the  total  stress  is  still  over  11,000  psi. 

From  point  "dw  to  point  "e*  the  sample  is  liquefied  and  the  total  stress 
equals  the  pore  pressure.  Also  note  that  the  hysteresis  in  the  total  stress 
curve  is  positive,  i.e.  energy  is  being  dissipated,  while  the  apparent 
hysteresis  in  the  pore  pressure  load-unload  cycle  is  negative.  The  negative 
hysteresis  results  from  the  rapid  drop  in  stress  on  the  solid  grains  as  the 
skeleton  unloads.  The  grains  tend  to  expand  rapidly  in  the  pore  fluid 
resulting  in  additional  pore  pressure  and  the  apparent  negative  hysteresis. 

At  point  MeM  the  pore  pressure  line  is  opened  and  shortly  thereafter 
total  stress  is  again  applied  and  the  sample  is  allowed  to  drain  and  con¬ 
solidate  until  the  effective  stress  reaches  its  pretest  value  at  point  "f" . 

The  strain  during  consolidation  to  the  original  effective  stress  is  nearly  80% 
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of  the  strain  reached  during  the  undrained  loading.  The  average  residual 
strain,  or  subsidence,  is  85%  of  the  peak  strain  reached  in  the  undrained 
loadings.  By  computing  the  peak  dynamic  strain  fields  from  explosive  detona¬ 
tions,  this  average  consolidation  can  then  be  used  to  estimate  the  portions  of 
total  crater  volumes  which  can  be  attributed  to  post-1 iquefaction  con¬ 
solidation. 

c.  Grain  Damage  as  a  Function  of  Stress  Path  and  Magnitude  -  Kim,  Blouin 
and  Timian,  1987  Section  7. 

A  systematic  study  of  grain  crushing  over  a  variety  of  stress  paths  to 
several  different  peek  mean  stresses  was  conducted.  Relating  the  measured 
grain  crushing  to  the  stress-strain  data  from  the  various  tests  provides 
insight  into  the  similarities  and  differences  batweon  the  different  types  of 
loadings  on  the  microscopic  levels,  and  particularly  into  the  differences  and 
similarities  between  the  microscopic  response  of  granular  materials  under 
shear  and  compressive  loadings.  Drained  hydrostatic  compression,  uniaxial 
strain,  and  triaxial  compression  tests  were  run  on  sieved  carbonate  beach  sand 
samples  having  a  uniform  grain  size  between  0.425  and  0.60  mm.  The  tests  were 
run  to  several  peak  mean  stresses,  thus  providing  a  variation  in  the  shear 
stress  component  of  the  loading  between  the  various  test  types. 

The  comparison  of  grain  damage  for  the  10,000  psi  mean  stress  tests, 
presented  in  Figure  1.2,  shows  a  moderate  to  high  degree  of  damage  with  36V  to 
49%  of  the  grains  broken  into  a  smaller  size  interval.  The  material  loaded  in 
hydrostatic  compression  sustained  the  least  damage,  with  36.4%  of  the  grains 
reduced  ih  size.  Damage  from  the  combined  hydrostatic  and  shear  loadings 
ranged  from  42.2%  for  the  TXC  test  to  49.1%  for  the  K0  test.  The  differences 
between  the  three  tests  is  in  proportion  to  the  shear  strain  energy  imparted 
to  each  sample,  with  the  hydrostatic  serving  as  the  zero  shear  baseline. 

In  order  to  more  clearly  define  the  role  of  shear  strain,  an  additional 
triaxial  compression  test  was  run  at  the  same  confining  pressure,  but  to  a 
much  higher  shear  strain  energy.  Severe  grain  damage  was  sustained  by  this 
triaxial  sample  loaded  under  the  same  8,400  psi  confining  pressure  ar  the 


10 


10.000  psi  mean  stress  sample,  but  carried  to  a  maximum  mean  stress  of 
15,350  psi  and  a  maximum  stress  difference  of  20,900  psi.  Nearly  80%  of  the 
grains  were  reduced  in  size,  with  over  half  the  total  post-test  weight  having 
grain  sizes  of  less  than  ,108  mm.  Shear  strain  energy  imparted  to  the  15,350 
psi  sample  was  about  an  order  of  magnitude  greater  than  that  of  the  10,000  psi 
triaxial  compression  test. 

In  order  to  concisely  and  quantitatively  describe  the  grain  damage  we 
developed  a  *echr.iquo  which  uses  the  pretest  and  post-test  log  mean  grain  sizes 
to  compute  a  grain  damage  factor,  Df.  The  grain  damage  factor  is  a  measure  of 
both  the  amount  and  severity  of  the  grain  damage.  A  sample  with  most  of  its 
grains  crushed  to  a  fine  powder  has  a  muc!,  larger  grain  damage  factor  than  a 
sample  with  the  same  percentage  of  its  grains  broken  into  only  slightly 
smaller  pieces.  As  shown  in  Table  1.2,  the  grain  damage  factor  ranges  from 
0.C11  for  the  1,000  psi  hydrostatic  sample  to  n.743  for  the  15,350  mean  stress 
triaxial  compression  sample.  The  maximum  possible  grain  damage  factor  is  1.0. 

1.2.2  Summary  of  Thworatlcal  Davmlopamnt.  code  Implementation  and  Numerical 
Analysis  From  Previous  Annual  Reports 

Most  of  the  theoretical  modal  development,  its  implementation  into  the 
codes,  and  verification  problems  and  parametric  analysis  using  the  codes  are 
summarized  in  Sections  2  through  5  of  this  report,  with  mu.h  of  the  detailed 
development  covered  in  Kim,  Blouin  and  Timian,  1986  and  1987.  The  theoretical 
formulations  for  the  numerical  code  NK0CP  for  the  prediction  of  undrained 
uniaxial  strain  and  hydrostatic  loadings  is  not  covered  in  these  sections. 

NK0CP  code  implementation  and  code  verification,  were  fui.y  described  by  Kim, 
Blouin  and  Timian,  1996  Section  5  and  Kim,  Slouin  and  Timien,  1987  Sections  2 
and  3  and  Appendices  A  and  8,  and  will  be  briefly  summarized  here. 

NKOCP  models  the  undrained  hydrostatic  and  uniaxial  strain  r*sponse  of 
saturated  porous  soils  and  rocks.  It  uses  the  measured  drained  skeleton 
properties  as  input  and  provides  the  total  stress  response,  effective  stress 
response  and  pore  pressure  response  as  output.  It  employs  an  incremental 
numerical  technique  to  model  the  nonlinear  fully  coupled  volumetric  response 
described  in  Section  2,  using  a  nonlinear  compressibility  of  fresh  water  or 
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sea  water  and  a  nonlinear  compressibility  model  for  the  solid  grains,  both 
described  by  Kim,  Blouin  and  Timian,  1987  Section  3. 

Comparison  of  an  undrained  uniaxial  strain  load-unload  cycle  on  saturated 
porous  limestone  computed  using  NKOCP  with  test  data  is  shown  in  Figure  1.3, 
parts  a  through  c.  The  agreement  between  the  numerical  calculation  and  the 
test  data  is  quite  good,  especially  considering  that  the  calculated  response 
was  developed  from  measured  drained  skeleton  response  on  a  similar,  but  not 
identical  rock  sample  and  from  constitutive  models  for  each  of  the  components 
which  were  developed  independently  of  the  test  data. 

During  the  current  year,  NKOCP  was  further  modified  to  incorporate  a 
nonlinear  unloading  capability  which  duplicates  the  actual  drained  skeleton 
unloading,  rather  than  the  original  bilinear  approximation  used  in  the  1987 
version.  This  new  capability  is  demonstrated  against  test  data  in  Section  7 
of  this  report. 

1.2.3  Summary  of  Experimental,  Theoretical  and  Numerical  Results  From 
Sections  2  Through  7 

a.  Section  2  -  Formulation  of  Field  Equations  for  the  Multiphase  Model 

In  Section  2  of  this  report,  the  theoretical  and  experimental  two  and 
three-phase  modeling  work  reported  in  the  two  previous  annual  reports  is  com¬ 
bined  with  work  performed  during  the  past  year  to  develop  the  final  for¬ 
mulation  of  the  field  equations  which  are  incorporated  into  MPDAP,  the 
multiphase  code  described  in  Sections  3  and  4. 

Features  and  advantages  of  the  new  formulations  over  the  previous 
multiphase  formulations  in  TPOAPII  include: 

1.  generalized  nonlinear  fluid  friction  equation  which  models  dyna¬ 
mic  pore  water  flow  In  both  the  laminar  and  turbulent  flow 
regimes;  and 

2.  fully  coupled  continuity  relationships  in  which  volumetric  strain 
compatibility  and  pressure/stress  equilibrium  between  the  pore 
fluid,  solid  grains  and  porous  skeleton  are  maintained. 
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The  six  governing  field  equations  used  in  MPDAP  are  summarized  at  the  end 
of  Section  2. 

b.  Section  3  -  Dynamic  Multiphase  Finite  Element  Formulations  for  MPDAP 

The  field  equations  in  Section  2  represent  the  fundamental  governing 
equations  for  an  infinitesimal  element  of  saturated  porous  medium.  Global 
equilibrium  equations  for  specified  boundary  conditions  are  derived  in  Section  3 
by  applying  the  principles  of  virtual  work  to  the  field  equations. 

Innovative  aspects  of  these  formulations  include: 

1.  formulation  of  the  global  equilibrium  equations  for  stresses  in 
the  bulk  medium  and  for  pore  fluid  flow  based  on  the  principals 
of  virtual  work  and  complimentary  virtual  work,  respectively;  and 

2.  use  of  pore  pressure  at  the  element  nodes  to  represent  relative 
motion  of  the  pore  fluid  which  eliminates  one  degree  of  freedom 
in  two-dimensional  calculations  and  two  degrees  of  freedom  in 
three-dimensional  calculations,  resulting  in  significant 
reductions  in  computational  running  time  and  storage  require¬ 
ments  . 

The  MPDAP  global  equilibrium  equations  are  summarized  in  Section  3.9. 

c.  Section  4  -  MPDAP  Verification  Problems 

The  features  and  capabilities  of  the  multiphase  code,  MPDAP,  are 
described  in  Section  3  and  a  number  of  verification  problems  presented  which 
compare  the  NPOAP  output  to  established  closed  form  solutions  of  simple 
problems  or  to  previous  TPOAP  solutions  for  more  complex  problems.  The  four 
verification  problems  presented  in  Section  4  include  comparison  of  MPDAP  solu¬ 
tions  with: 

1.  the  closed  form  solution  for  an  undrained  uniaxial  strain 
loading; 

2.  the  closed  form  solution  for  spherical  elastic  wave  propagation 
in  a  single  phase  medium; 
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3.  Terzaghi's  closed  form  solution  for  quasi-static  flow  and  con¬ 
solidation  under  a  uniaxial  strain  loading;  and 

4.  TPDAPII  solutions  for  one-dimensional  wave  propagation  resulting 
from  simulated  explosive  loadings  of  both  a  saturated  soil  and 
rock. 

d.  Section  5  -  Numerical  and  Theoretical  Treatment  of  Waves  of  the 

First  and  Second  Kind 

In  Section  5,  theoretical  development  and  coding  of  a  solution  for  wave 
propagation  in  saturated  porous  elastic  media  is  reviewed.  This  work, 
reported  in  detail  by  Kim,  Blouin  and  Timian,  1967  Sections  4  and  5,  resulted 
in  the  code  TWAVE  which  uses  the  closed  form  solutions  to  rapidly  compute 
wavespeeds  and  damping  over  a  wide  range  of  material  properties.  Significant 
findings  from  TWAVE  include  the  following: 

1.  Existence  of  two  types  of  compressions!  waves  in  saturated  porous 
media;  conventional  compress ional  waves  analogous  to  those  in 

a  single  phase  material,  termed  waves  of  the  first  kind,  and 
much  slower  highly  damped  compressional  waves  termed  waves  of  the 
second  kind. 

2.  Wavespeeds  of  both  types  of  waves  have  a  lower  and  upper  bound 
wavespeed  which  is  dependent  on  the  excitation  frequency  and/or 
the  material  permeability.  Wavespeed  increases  with  increasing 
excitation  frequency  and  with  increasing  permeability. 

3.  In  waves  of  the  first  kind,  the  rate  of  wavespeed  increase  is 
greatest  at  frequencies  and  permeabilities  around  which  damping 
is  a  maximum. 

In  Section  5  of  this  report,  TPDAPII  solutions  for  wave  propagation 
having  skele.Jn  properties  and  permeabilities  chosen  to  highlight  the 
wavespeed  and  damping  differences  are  compared  to  the  TWAVE  solutions  and 
analyzed  to  determine  the  factors  governing  the  damping  variations  and 
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wavespeed  changes  and  to  determine  the  phenomena  behind  propagation  of  waves 
of  the  second  kind.  Conclusions  from  this  comparison  include  the  following: 

1.  Both  the  material  skeleton  and  pore  water  are  in  compression 
during  passage  of  waves  of  the  first  kind.  Pore  fluid  friction 
is  at  a  maximum  where  damping  and  wavespeed  increases  are  at  a 
maximum.  At  the  lower  bound  wavespeed,  there  is  no  relative 
motion  between  the  pore  water  and  the  porous  skeleton. 

2.  Waves  of  the  second  kind  appear  to  be  associated  with  a  surge  of 
pore  fluid  moving  through  the  skeleton.  Pressure  in  the  pore 
water  is  compressive,  while  stress  in  the  porous  skeleton  is  ten¬ 
sile  during  passage  of  waves  of  the  second  kind.  The  pore  water 
is  moving  in  the  direction  of  wave  propagation  while  the  porous 
skeleton  is  moving  in  the  opposite  direction. 

e.  Section  6  -  Experimental  Evaluation  of  Fluid  Flow  in  Ducts  and  Soils 

A  series  of  fluid  flow  test  data  is  reported  by  Kim,  Blouin  and  Timian, 
1987  Section  8  and  in  Section  6  of  this  report.  Flow  velocities  of  up  to  3700 
in/s,  fluid  acceleration  to  nearly  1000  Qs  and  pressure  gradients  of  up  to  350 
psi/in  were  achieved  in  some  of  the  dynamic  flow  tests.  Results  of  this  study 
were  as  follows: 

1.  Recommendation  of  a  fluid  flow  equation  which  models  flow  in 
both  the  laminar  and  turbulent  flow  regimes; 

2.  Development  of  an  apparatus  and  experimental  techniques  for  the 
determination  of  the  flow  coefficients  in  the  above  equations 
for  porous  soils  and  rocks;  and 

3.  Demonstration  that  pore  fluid  flow  for  explosive  loadings  will  be 
largely  in  the  turbulent  flow  regime  where  the  pore  pressure 
gradient  is  proportional  to  the  square  of  the  pore  fluid  velocity 
relative  to  the  skeleton.  It  was  also  demonstrated  that  Biot's 
theoretical  formulation  for  the  increase  in  fluid  friction  with 
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increasing  excitation  frequency  can  be  ignored  in  the  turbulent 
regime. 

f.  Section  7  -  Two-Phase  Response  of  Saturated  Limestone 

In  Section  7,  a  series  of  drained  and  undrained  test  data  on  Indiana 
limestone  are  analyzed  using  NKQCP  and  two-phase  constitutive  properties. 
Orained  and  undrained  hydrostatic,  uniaxial  strain,  triaxial  compression  under 
constant  confining  pressure  and  specified  strain  path  tests  were  conducted. 
Analysis  of  the  four  test  types  showed  that: 

1.  Undrained  hydrostatic  and  uniaxial  strain  response  was  in 
excellent  agreement  with  calculations  using  the  revised  version 
of  NKOCP.  Total  stress,  effective  stress  and  pore  pressure 
response  closely  matched  the  test  data  over  the  entire 
load-unload  cycle. 

2.  Undrained  triaxial  strength  and  deformation  prior  to  failure 
were  in  good  agreement  with  predictions  based  on  effective 
stress  theory  and  two-phase  models.  However,  strengths  and 
pore  pressure  response  during  the  latter  stages  of  shearing 
appear  to  contradict  the  response  predicted  by  effective  stress 
theory  and  two-phase  models.  We  believe  this  apparent  contra¬ 
diction  is  due  to  inherent  shortcomings  of  the  test  itself, 
caused  by  lateral  constraint  from  the  steel  end  caps.  Use  of 
such  test  data  could  lead  to  gross  inaccuracies  in  modeling  in 
situ  response. 

3.  Drained  and  undrained  strain  path  tests  that  imposed  a  variety 
of  strain  paths  on  the  test  specimens  following  uniaxial  strain 
compression  will  prove  to  be  a  real  challenge  to  material 
modelers.  These  strain  path  tests,  which  are  typical  of  strain 
paths  from  explosive  loadings,  point  up  the  need  for  an  advanced 
single  element  two-phase  code  which  can  exercise  various 
saturated  skeleton  models  over  arbitrary  strain  and  stress 
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paths.  Test  results  were  generally  consistent  with  response  that 
would  be  predicted  with  two-phase  models,  though  the  effective 
stress  response  may  be  difficult  to  match  with  existing  skeleton 
models.  The  test  results  highlight  the  dependence  of  strength 
on  the  strain  path  (as  predicted  by  effective  stress  theory  and 
two-phase  models)  and  clearly  demonstrate  the  inadequacy  of  using 
equivalent  single-phase  calculations  for  predicting  response 
to  explosive  loadings  in  saturated  porous  materials. 


17 


Table  1.1  Bulk  modulus  values  obtained  from  grain  compressibility  tests 


Material 

Bulk 

Modulus 

(psi  x  106) 

(MPa  x  104) 

Steel  Ball  Bearings 

22.2 

15.31 

Quartz  Sand 

5.2 

3.58 

Solid  Limestone 

11.0 

7.59 

Eriewetak  Beach  Sand 

9.1 

6.27 

Silt-Sand-Gravel  from  KAM-2 

5.6 

3.86 

Ground  Silt-Sand-Gravel  from  KAM-2 

6.1 

4.21 

Vugular  Limestone 

5.7 

3.93 

Cemented  Material  from  XSA-2 

6.1 

4.21 

Table  1.2a.  Summary  of  post-test  grain  size  distribution  and  grain  damage;  uniform 
Enewetak  beach  sand,  Gs  =  2.78  qm/cwfi. 
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Grain  Size  (mm) 

Figure  1.2.  Comparison  of  post-test  grain  size  and  weight  fraction  distribution 
from  samples  loaded  to  10,000  psi  and  15,350  psi  mean  stress. 
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numerical  simulation. 
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Figure  1.3c.  Comparison  of  effective  axial  stress  response  of  saturated  porous 
limestone  with  numerical  simulation. 


SECTION  2 


FORMULATION  OF  FIELD  EQUATIONS 
FOR  THE  MULTIPHASE  MODEL 


2.1  INTRODUCTION 

The  initial  theoretical  formulations  to  be  incorporated  into  the 
multiphase  code  MPDAP  were  originally  presented  by  Kim,  Blouin  and  Timian, 

1986  and  1987.  During  the  past  year  the  original  formulations  have  been 
revised  to  provide  more  efficient  and  accurate  computational  algorithms 
including  the  modification  of  the  dynamic  flow  equation  based  on  the  results 
of  our  dynamic  flow  tests  described  in  Section  6. 

The  revised  theoretical  formulations  have  been  incorporated  into  the 
MPOAP  code  during  the  past  year  and  verification  problems  have  been  run. 
Features  and  advantages  of  the  new  formulations  over  previous  multiphase  for¬ 
mulations  include  the  following: 

1.  Generalized  nonlinear  fluid  friction  equation  which  models  dynamic 
pore  water  flow  in  both  the  laminar  and  turbulent  flow  regimes; 

2.  Fully  coupled  continuity  relationships  in  which  volumetric  strain 
compatibility  and  pressure/stress  equilibrium  between  the  pore  fluid, 
solid  grains  and  porous  skeleton  are  maintained; 

3.  Formulation  of  the  global  equilibrium  equations  for  stresses  in  the 
bulk  medium  and  for  pore  fluid  flow  based  on  the  principals  of 
virtual  work  and  complimentary  virtual  work,  respectively;  and 

4.  Use  of  pore  pressure  at  the  element  nodes  to  represent  relative 
motion  of  the  pore  fluid  which  eliminates  one  degree  of  freedom  in 
two-dimensional  calculations  and  two  degrees  of  freedom  in  three- 
dimensional  calculations,  resulting  in  significant  reduction'  in 
computational  running  time  and  storage  requirements. 
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2.2  THEORETICAL  FORMULATIONS 


2.2.1  Notation 

Note  that  positive  signs  have  been  used  for  elongation  and  tension.  A 
comma  denotes  differentiation  with  respect  to  the  subsequent  indices  and  the 
superposed  dot  denotes  time  rate.  Prime  indicates  effective  stress  or 
pressure. 

<u)  :  skeleton  displacement 

{U}  :  absolute  fluid  displacement 

<w)  :  apparent  fluid  displacement  relative  to  the  sol io  skeleton 

(or)  :  total  stress 

(o')  :  effective  stress 

p  :  total  pressure 

p'  :  effective  pressure 

n  :  pore  fluid  pressure 

:  pore  fluid  pressure  gradient  vector 

{ e>  :  skeleton  strain 

ev  :  skeleton  volumetric  strain 

€f  :  pore  fluid  volumetric  strain 
eg  ;  solid  grain  volumetric  strain 
ep  :  volumetric  diffusion  of  pore  fluid 
{u)e  ;  element  nodal  skeleton  displacement  vector 
{ rr) e  :  element  nodal  pore  fluid  pressure 

(ul  :  global  nodal  skeleton  displacement  vector 
in)  :  global  nodal  fluid  pressure 

<T>  :  applied  boundary  traction 

0  :  specified  boundary  flow  velocity  (flux) 

(b)  :  body  force  vector  (generally  equals  gravity  force) 

k  :  Darcy's  coefficient  of  permeability 

[DeP]  :  elasto-plastic  stress-strain  matrix  for  skeleton 

(1)  ;  unit  vector  <1)T  =  <1  l  l  Q  00> 

n  ;  porosity 

Cf  :  pore  fluid  compressibility 
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caw 


Km 

Ks 

Kf 

Ms 

P 

Pd 

Pf 

p9 

Vt 

Yf 


Of 

e 

[Mt] 

[KT1 

[C] 

[E] 

[H] 

<F> 

<R> 

<Q> 


:  compressibility  of  air-water  mixture 
:  compressibility  of  solid  grains 

:  compressibility  of  soil-water  mixture  with  zero  effective  stress 
:  bulk  modulus  of  soil  water  mixture  with  zero  effective  stress 
:  bulk  modulus  of  skeleton 
:  bulk  modulus  of  pore  fluid 
:  constrained  modulus  of  skeleton 
:  bulk  mass  density  of  mixture 
:  dry  density  of  skeleton 
;  fluid  mass  density 
:  mass  density  of  solid  grains 
:  mass  of  solid  grains  in  volume  V^. 

:  total  skeleton  volume 
:  unit  weight  of  the  pore  fluid 

:  mass  increment  factor  (approximation  of  Biot's  fluid  viscosity 
parameter) 

s  Ward's  fluid  friction  coefficient  for  turbulent  flow 
:  parameters  in  Newmark's  fi  numerical  time  integration  method 
:  parameter  in  Wilson's  9  numerical  time  integration  method 
:  Kronecker's  delta 
:  Mass  Matrix 

:  tangent  skeleton  stiffness  matrix 

:  coupling  matrix  between  solid  skeleton  and  pore  fluid 

:  pore  fluid  compressibility  matrix 

:  fluid  friction  energy  dissipation  matrix 

s  nodal  force  vector 

:  internal  resistance  force  vector 

:  equivalent  boundary  flow  vector 


2.2.2  Field  Equations 


Effective  Stress  Law 

Terzaghi's  effective  stress  equation  is  fundamental  to  the  development  of 
the  fully  coupled  model.  It  relates  the  total  applied  stress,  a,  to  the  pore 


31 


pressure,  n,  and  the  effective  stress,  o',  according  to 


^i  j  =  0,1  ij  +  5ij7r 


(2-1) 


where  CTij  =  total  stress 

o'  .j  j  =  effective  stress 
5jj  =  Kronecker's  delta 
fiij  =  0  if  i  *  j 
$jj  =1  if  i  =  j 

Constitutive  Equation  for  Skeleton  Deformation 


The  deformation  of  the  porous  skeleton  is  related  to  the  applied  effec¬ 
tive  stress  and  the  pore  pressure  acting  on  the  solid  grains.  The  stress- 
strain  relationship  is  given  by 


{da'}  =  [D«P]  ({de}  -  23  |l}dn) 

3 


(2-2) 


The  last  term  in  Equation  2-2  is  the  strain  in  the  skeleton  resulting  from 
compression  of  the  solid  grains  by  the  pore  pressure. 

Continuity  Equation  of  Pore  Fluid  Flow 

The  continuity  equation  for  pore  fluid  flow  is  derived  from  mass  conser¬ 
vation  relationships.  The  volumetric  strain  of  the  pore  fluid,  Cf,  is  given 

by 

dfif  a  -  f|£f  a  Cf  dn  (2-3) 

pf 

where  Cf  =  pore  fluid  compressibility 

7T  =  pore  fluid  pressure 
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(2-4) 


where  Cg  =  bulk  compressibility  of  solid  grains 

p'  =  effective  mean  pressure 

The  dry  density,  pd,  is  given  by 

pd  =  =  (1  -  n)pg  (2-5) 

where  mg  is  the  mass  of  the  solid  grains  in  skeleton  volume  V*.  The  change  in 
dry  density  is  given  by 


dpd  3  -Pd  dev  (2-6) 

where  ev  =  is  the  volumetric  strain  of  the  skeleton.  Differentiating 
Equation  2-5  with  respect  to  n  and  pg  gives 

dpd  ■  (1  -  n)  dpg  -  Pgdn  (2-7) 


Equating  2-6  and  2-7  yields 


(2-8) 


Conservation  of  mass  for  the  pore  fluid  within  a  specified  initial  volume 
of  saturated  porous  material  is  given  by 

n  pf  Vt  =  n  Pf  Vt  (2-9) 

where  as  illustrated  in  Figure  2.1,  the  terms  to  the  left  of  the  equal  sign 
represent  the  fluid  mass  under  the  initial  conditions  and  the  terms  to  the 
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right  represent  the  same  fluid  mass  under  deformed  conditions.  Equation  2-9 
may  be  expressed  in  infinitesimal  incremental  form  as 

npfVt  =  (n  +  dn)(pf  +  dpp)(l  +  dep)  Vt  (2-10) 


where 


£p  =  volumetric  diffusion  of  pore  fluid  as  depicted  in  Figure  2.1. 
Solving  Equation  2-10  for  dep  and  discarding  second  order  terras  yields 


(2-11) 


Equation  2-11  is  combined  with  Equation -.2-8  by  elimination  of  dn  to  yield 


(1  -  n)dev  ♦  ndep  ♦  (1  -  n)  ♦  n  «  0 


(2*12) 


Combining  Equations  2-3  and  2-4  with  2-12  gives 


n(dep  -  dev)  +  dev  -  ^  dn  -  Cgdp*  *»  0  (2-13) 


where  Kra  is  the  bulk  modulus  of  the  solid/fluid  mixture  which  is  expressed  by 


=  nCf  ♦  (1  -  n)Cg 

The  change  in  effective  mean  pressure  is  given  by 


(2-14) 
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dp'  =  Ks  (dev  -  Cgdn) 


(2-15) 


Substituting  Equation  2-15  into  2-13  gives 


The  differential  equation  of  notion  governing  the  bulk  mixture  is 
expressed  by  equating  the  stress  gradient  to  the  inertial  resistance  as 


°ij.j  *  (*  ’  ">Ps  ui  *  °Pf  ^i  (2-H) 

o^j#j  is  the  total  stress  gradient  applied  to  an  infinitesimal  element 
of  saturated  material  at  some  given  tine,  o'ijj  is  expressed  in  tensor 
notation  and  represents  the  stress  gradient  in  each  of  three  mutually  perpen¬ 
dicular  coordinates  (e.g,  see  Handle son,  1968).  For  instance,  in  the  x 
direction, 

do  da  3a 

axj.j  8  *  “3^  3  U  -  n)psux  ♦  nPf  (*x  (2-18) 

The  term  (1  -  n)ps  is  the  mass  of  the  soil  skeleton  per  unit  volume  of 
saturated  material,  where  n  is  the  porosity  and  py  is  the  mass  density  of  the 
solid  grains,  u^  is  the  displacement  of  the  skeleton  in  the  i  direction  and 
ii,  is  the  acceleration  of  the  skeleton  in  the  i  direction.  The  term  npf  is  the 
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mass  of  pore  fluid  per  unit  volume  of  saturated  material  where  pf  is  the  mass 
density  of  the  pore  fluid.  is  the  absolute  displacement  of  the  pore  fluid 
in  the  i  direction. 

The  bulk  mass  density  of  the  saturated  material,  p,  is  given  by 


p  =  (1  -  n)ps  +  npf 


(2-19) 


Substitution  of  the  value  for  (1  -  n)ps  from  Equation  2-19  into  Equation  2-17 
gives 


ffij,j  B  "  "PfK'i  *  "Pf^i 


(2-20) 


A  tern  is  introduced  which  is  the  apparent  fluid  displacement  in  the  i 
direction  relative  to  the  soil  skeleton, and  is  given  by 

w*  «  n(Ui  “  Uf)  (2-21) 

In  seepage  problems,  w^,  is  referred  to  as  the  discharge  displacement.  It 
describes  the  discharge  of  fluid  through  a  soil  mass  of  unit  area.  The 
discharge  velocity,  or  apparent  relative  velocity,  Wj,  between  the  soil  par¬ 
ticles  and  pore  water  is  the  velocity  of  water  in  a  discharge  duct  of  unit 
area  needed  to  maintain  the  actual  relative  velocity  in  the  porous  soil  of  the 
same  unit  area.  The  actual  relative  velocity  between  the  skeleton  and  the 
pore  water  is  given  by  w^/n.  Finally,  is  the  apparent  relative  accelera¬ 
tion  between  the  soil  skeleton  and  pore  water  given  by 

w^  *  n(Ui  -  Uj)  (2-22) 

Equation  2-20  can  be  expressed  in  terms  of  the  apparent  relative  fluid  acce¬ 
leration  as  simply 
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Equation  of  Motion  for  Pore  Fluid 


Over  the  past  several  years  our  equation  describing  the  motion  of  the 
pore  fluid  relative  to  the  skeleton  has  undergone  several  evolutionary  revi¬ 
sions.  The  original  version,  described  by  Kim  and  Blouin,  1984,  utilized  an 
approximation  to  Biot's  fluid  friction  equations  given  by 


ff,i  .  pfUj  +  +  ir  "i  +  «i 


Inertial 

Components 


Frictional 

Components 


(2-24) 


where  n ^  =  pore  pressure  gradient 

*  absolute  skeleton  acceleration 
w^  »*  apparent  fluid  velocity  relative  to  the 
skeleton 

=  apparent  fluid  acceleration  relative  to 
the  skeleton 

pf  =  pore  fluid  mass  density 
yf  ~  unit  weight  of  the  pore  fluid 
n  =  porosity 

k  =  Darcy's  coefficient  of  permeability 
r  =  mass  increment  factor 

The  first  two  terms  to  the  right  of  the  equal  sign,  the  inertial  com¬ 
ponents,  represent  the  portion  of  the  pressure  gradient  resulting  from  acce¬ 
leration  or  deceleration  of  the  pore  fluid.  The  last  two  terms  represent  the 
portion  of  the  pore  pressure  gradient  due  to  fluid  friction  associated  with 
the  relative  motion  between  the  pore  fluid  and  solid  skeleton.  The  two  fluid 
friction  terras  are  based  on  Biot's  work  and  are  proportional  to  the  relative 
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fluid  velocity  and  acceleration,  respectively.  The  first  friction  terra  is 
identical  tc  Darcy's  law  for  steady  state  flow  conditions.  The  second  fric¬ 
tion  term  is  a  general ization  of  Biot's  frequency  dependent  friction  term. 

Biot  expressed  the  pore  pressure  gradient,  as 

JT,j  =  PfU-j  +  D-j  (2-25) 


where  PfU^  is  the  inertial  force  per  unit  volume  of  pore  fluid  and  D-j  repre¬ 
sents  the  viscous  friction  force  between  the  pore  fluid  and  the  soil  skeleton 
per  unit  volume  of  pore  fluid.  Solving  Equation  2-22  for  U-j  and  substitution 
into  Equation  2-25  gives 


JT#i  a  El  Wj  *  PfU-j  *  D-j 


(2-26) 


Biot  showed  that  the  viscous  friction  term,  D-j ,  is  a  function  of  the 
excitation  frequency,  w,  the  pore  geometry,  the  dynamic  viscosity,  p,  and  the 
apparent  relative  velocity  between  the  pore  fluid  and  the  skeleton,  Wj.  In  an 
actual  soil  the  flow  of  pore  water  would  follow  very  complicated  paths  which 
are  difficult  to  describe.  These  flow  paths  would  involve  numerous  variations 
in  direction  and  in  cross  sectional  area.  Biot  employed  models  of  the  flow 
paths  which  are  gross  simplifications  of  the  actual  paths.  He  assumed  two 
simple  flow  geometries ?  flow  through  a  series  of  parallel  circular  ducts  and 
flow  through  a  series  of  parallel  flat  ducts.  Flow  conditions  in  the  flat 
duct  are  depicted  schematically  in  Figure  2.2. 


For  dynamic  laminar  flow,  Biot  (1956)  derived  the  exact  expression  for 
the  viscous  friction  term,  D-j,  for  both  circular  and  flat  ducts  which  is  given 
by 

»i  =  Yrj^1  *i  (2-27) 


where  yf  is  the  unit  weight  of  pore  fluid,  F(x)  is  the  viscous  friction 
correction  factor,  and  k  is  the  coefficient  of  permeability.  F(k)  is  a 
complex  function  given  by 


38 


F(tc)  =  fX(K)  +  i  f2(ic) 


(2-28) 


The  magnitudes  of  the  real  part,  f ^ (k) ,  and  imaginary  part,  f 2 (K) ,  are  plotted 
in  Figure  2.3  as  a  function  of  nondimensional  parameter,  k,  which  is  defined 
as 


K  a  P  JC0 


(2-29) 


where 


*o 


=  (JL) 

vn  g' 


w  k 


(2-30) 


In  Equations  2-29  and  2-30,  g  is  the  gravitational  acceleration,  w  is  the 
excitation  frequency,  and  the  factor,  r,  is  the  constant  which  is  dependent  on 
the  shape  of  the  flow  path.  For  Generalized  Darcy's  Flow  (also  called 
Poiseuille  Flow),  f ^ (ic)  and  f 2(k)  are  independent  of  tc  and  r  =  0.  For  the 
circular  duct,  r  is  unity.  And  for  the  flat  duct,  the  value  of  r  is  approxi¬ 
mately  equal  to  /i73. 

An  approximation  of  Equation  2-27  was  developed  by  Kim  and  Blouin  (1984) 
and  is  given  by 


Di 


3 


Vf  ^ 

IT  * 


f  E!  r  iii 
n  1 


(2-31) 


where  tf  »  unit  weight  of  pore  fluid 

r  0  empirical  mass  increment  factor  which  treats 
the  dynamic  fluid  friction  as  though  it 
was  an  increase  in  inertial  mass  of  r 

The  theoretical  value  of  r  is  1/3  for  the  circular  duct  and  1/5  for  the  flat 
duct.  Equation  2-31  is  a  good  approximation  when  the  nondimensional  para¬ 
meter,  k,  is  less  than  2. 


Substitution  of  Equation  2-31  into  Equation  2-26  gives 


nti  ■  £f  (1  +  r)  Wi  +  pfUi 


*P 


W4 


(2-32) 


which  has  been  rewritten  as  Equation  2-24  to  separate  the  frictional  terms 
from  inertial  terms. 

In  order  to  better  understand  the  influence  of  Biot's  dynamic  frictional 
resistance  the  normalized  relative  flow  velocity  distribution  in  a  flat  duct 
given  by  Biot  (1956)  is  plotted  in  Figure  2.4  as  a  function  of  the  fluid 
viscosity  and  excitation  frequency.  The  shape  of  the  velocity  profile  is  a 
function  of  the  nondimensional  parameter  0  which  is  defined  as 


where  aj  =  half  height  of  flat  duct 

w  =  excitation  frequency 
v  =  kinematic  viscosity 

Figure  2.4  illustrates  the  influence  of  excitation  frequency  on  the  shape 
of  the  flow  velocity  distribution  in  the  flat  duct  for  a  constant  fluid  visco¬ 
sity.  For  low  excitation  frequencies  the  distribution  is  parabolic,  the  same  as 
that  for  steady  state  laminar  flow.  As  the  excitation  frequency  increases, 
the  velocity  distribution  is  pinched  toward  the  center  of  the  duct,  indicating 
widening  static  boundary  layers  along  the  duct  walls. 

The  next  addition  to  the  equation  of  motion  for  the  pore  fluid  was  an 
additional  term  to  account  for  frictional  energy  losses  during  turbulent  flow 
conditions.  Figure  2.5  shows  theoretical  velocity  distributions  in  a  flat 
duct  in  both  the  laminar  and  turbulent  flow  regimes.  Laminar  flow  occurs  at 
Reynolds  numbers  of  less  than  2000.  For  steady  state  laminar  flow  the  velocity 


40 


u 

u 


i 


» 


distribution  is  parabolic  about  the  center  of  the  duct  (Pouiseuille  flow).  For 
nonsteady  state  flow,  Biot's  theoretical  analysis  indicates  a  change  in  the 
flow  velocity  profile  such  as  sown  in  Figure  2.4  and  2.5.  In  this  flow  state 
fluid  friction  increases,  but  flow  remains  laminar. 

When  the  Reynolds  number  is  greater  than  approximately  2000,  the  flow 
becomes  turbulent  with  a  velocity  distribution  in  the  direction  of  flow  simi¬ 
lar  to  that  shown  in  Figure  2.5,  where  there  is  a  very  sharp  velocity  gradient 
adjacent  to  the  duct  walls.  The  flow  lines  in  the  turbulent  regime  are  no 
longer  parallel  and  there  are  random  particle  motions  transverse  to  the  direc¬ 
tion  of  flow. 

As  described  in  Section  6,  we  have  experimentally  verified  a  generalized 
turbulent  flow  relationship  of  the  form 


=  aw-j  +  bjb-j2  +  cw.(  + 


(2-34) 


where  it,)  »  pore  pressure  gradient 

w^  =  apparent  flow  velocity  relative  to  the 
solid  skeleton 

w^  =  apparent  relative  pore  fluid  acceleration 
a,b  &  c  =  flow  parameters  which  are  a  function  of  microscale 
parameters  including  pore  sizes  and  geometries, 
overall  porosity,  viscosity,  etc. 

Pf  =  fluid  mass  density 
ii<j  a  skeleton  acceleration 


41 


In  the  laminar  flow  regime  the  parameters  in  Equation  2-34  are  given  by 


b  »  0  \  (2-35) 

c  =  (1+r)  ££ 

where  y^  =  fluid  unit  weight 

Pf  =  fluid  mass  density 
k  =  coefficient  of  permeability 
r  »  mass  increment  factor  which  approximates 
the  influence  of  excitation  frequency 
(r  =  1/5  for  flat  ducts  and 
r  =  1/3  for  circular  ducts) 

In  the  turbulent  flow  regime  the  flow  parameters  of  Equation  2-34  are 
given  by 


where 


(2-36) 


c  arr  ’ 

a  Ward's  turbulent  flow  constant  which  is 
a  function  of  the  pore  geometry  and  pore  fluid 
viscosity  and  which  is  determined  experiment¬ 
ally  from  steady  state  turbulent  flow  tests 
(see  Kira,  Blouin  and  Tinian,  1987). 


Note  that  in  the  tu.  bulent  regime  the  random  transverse  particle  motions 
result  in  elimination  of  the  mass  increment  factor  r  from  the  flow  parameter 


42 


c.  Section  6  presents  validation  of  Equation  2-36  for  flow  in  a  flat  duct 
with  velocities  up  to  4000  in/s  and  fluid  accelerations  to  nearly  1000  g's. 
Flow  parameters  a  and  b  have  also  been  validated  for  steady  state  turbulent 
flow  in  porous  soils. 

The  experimental  results  in  Section  6,  along  with  the  previous  experimen¬ 
tal  work  (Kim,  Blouin  and  Timian,  1986  and  1987),  demonstrate  that  the  tur¬ 
bulent  pore  water  flow  regime  will  dominate  the  fluid  friction  response  in' 
explosive  loadings  of  saturated  soils  and  rocks  at  stress  levels  of  interest 
to  civil  engineering  problems.  Turbulent  flow  develops  nearly  simultaneously 
with  arrival  of  the  dynamic  stress  wave  from  an  adjacent  explosion.  Laminar 
(or  Oarcy)  flow  is  relevant  only  at  late  times  after  passage  of  the  dynamic 
stress  waves.  The  experimental  and  theoretical  results  demonstrate  that 
Biot's  theoretical  work  (1956,  1962A,  1962B),  which  is  used  by  some  investiga¬ 
tors  to  describe  transient  flow,  is  not  applicable  in  the  turbulent  flow 
regime. 

Partial  Saturation 

In  numerous  geologic  settings  the  soil  or  rock  is  not  fully  saturated. 
Rischbieter  et  al.  (1977)  demonstrated  that  even  a  minute  amount  of  entrapped 
air  drastically  alters  the  pore  pressure  response  in  multiphase  porous 
materials.  Thus,  a  complete  treatment  of  multiphase  media  should  include  the 
capability  of  calculating  stress  wave  propagation  and  pore  fluid  response  in 
three-phase  porous  materials. 

Kim  (1982)  developed  a  unique  formulation  for  the  compressibility  of  the 
air-water  mixture  in  partially  saturated  porous  media.  This  formulation  has 
been  extensively  applied  in  quasi-static  problems  and  verified  against  experi¬ 
mental  data. 

The  compressibility,  Caw,  of  tha  air-water  mixture  in  partially  saturated 
media  is  given  by 
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(2-37) 


Caw  =  “  Sq  +  HCS0) 


^ao 

(n  +  T)2 


where  S0  =  initial  degree  of  saturation 

Hc  =  coefficient  of  solubility  (Henry's  constant) 
rrao  s  initial  pore  air  pressure  (absolute) 
n  =  current  pore  water  pressure  (absolute) 

T  =  pressure  difference  between  the  air  and  pore 
water  due  to  surface  tension 

As  indicated  by  Equation  2-37,  the  compressibility  of  the  air-water  mixture  is 
nonlinear  with  respect  to  the  current  pore  water  pressure.  This  relationship 
has  been  shown  to  be  applicable  when  the  degree  of  pore  water  saturation  is 
above  approximately  85%,  and  the  air  water  mixture  is  thought  to  exist  in  an 
occluded  state  with  the  air  contained  as  small  bubbles  within  the  pore  fluid. 

The  pressure  difference,  T,  in  Equation  2-37  is  determined  experimentally 
from  undrained  hydrostatic  tests  on  partially  saturated  samples  with  a  known 
degree  of  saturation  (Kim,  1982). 

The  treatment  of  partially  saturated  materials  utilizes  the  same  flow 
field  equations  as  described  in  the  preceding  discussion.  However,  the 
compressibility  of  the  pore  fluid,  Cf,  is  given  by  Equation  2-37.  In  other 
words,  the  porous  skeleton  is  assumed  to  be  filled  with  an  equivalent  fluid 
which  has  the  same  compressibility  as  the  air-water  mixture. 

At  some  value  of  pore  fluid  pressure  the  air  bubbles  undergo  a  complete 
collapse  and  the  compressibility  of  the  air-water  mixture  becomes  approxima¬ 
tely  equal  to  the  compressibility  of  the  pore  water. 

Summary  of  Field  Equations 

A  summary  of  all  the  field  equations  implemented  in  MPDAP  includes: 
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1.  Effective  Stress  Law, 


2.  Skeleton  Deformation  Relationship, 


3.  Pore  Fluid  Continuity  Equation, 


where  a  = 

and  b  «  0 

c  =  (1+r)  £f 

for  laminar  flow  conditions,  and 
b  = 

‘■P 

for  turbulent  flow  conditions 


(2-1) 


(2-2) 


(2-16) 


(2-20) 


(2-34) 
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Conservation  of  Fluid  Mass 


n  pf  Vt  =  n*  pf'  Vt' 


Vt  =  apparent  fluid  volume  before  compression 

Vt'  =  (1  +  ep)  Vt:  apparent  fluid  volume  after  compression 

ev  =  volumetric  strain  of  porous  skeleton 

6p  =  volumetric  diffusion  of  pore  fluid 


Before  Compression  After  Compression 


Figure  2.1.  Schematic  illustration  of  conservation  of  pore 
fluid  mass  in  saturated  porous  materials. 


Figure  2.2  Schematic  view  of  fluid  flow  in  a  flat  duct. 
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Figure  2.5.  Velocity  distributions  in  a  flat  duct. 


SECTION  3 


DYNAMIC  MULTIPHASE  FINITE  ELEMENT  FORMULATIONS  FOR  MPDAP 

3.1  INTRODUCTION 

The  field  equations  in  Section  2  represent  the  fundamental  governing 
equations  for  an  infinitesimal  element  of  saturated  porous  medium.  Using  the 
finite  element  method,  the  values  of  the  field  variables  at  any  point  within 
the  element  can  be  computed  from  the  element  nodal  values.  Global  equilibrium 
equations  for  specified  boundary  conditions  can  then  be  derived  by  applying 
the  principles  of  virtual  work.  These  global  equations  are  incorporated  into 
MPOAP  as  described  in  this  section. 

As  noted  in  the  previous  section,  pore  pressure  at  the  element  nodes  is 
used  to  compute  pore  pressure  gradients,  which  are  in  turn  used  to  compute  the 
relative  flow  of  the  pore  fluid.  This  approach  is  used  in  place  of  the  more 
conventional  technique  whereby  relative  fluid  displacement  vectors  are  defined 
at  each  element  node.  This  eliminates  one  degree  of  freedom  per  node  in  two- 
dimensional  calculations  and  two  degrees  of  freedom  per  node  in  three- 
dimensional  calculations,  resulting  in  large  reductions  in  running  time  and 
storage  requirements. 

3.2  SPATIAL  DISCRETIZATION 

Within  each  element,  field  variables  can  be  expressed  in  terms  of  the 
element  nodal  values  using  the  assumed  shape  functions. 

M  =  IN]  {0|e 
U)  =  {B)  |0je 

(3-1) 

n  =  <G>  {«}* 

(n,i)  *  [a]  (Jl)e 

where  the  above  notation  is  described  in  Section  2.2.1. 
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3.3  INCREMENTAL  FORM  OF  TOTAL  STRESS  VECTOR 


The  total  stress  vector  at  time  step  n  can  be  expressed  as 

W  =  { ^n- 1 }  +  I Acr  1  n  1  *  U}  (3-2) 

Substitution  of  Equation  2-2  into  Equation  3-2  yields 


W  ■  (oh-i)  +  [o^](*««}  *  (li)  -  [o*p]|i|)to„ 


(3-3) 


3.4  APPROXIMATION  OF  ABSOLUTE  FLUID  MOTION 

In  the  previous  section,  we  presented  Equation  2-35  describing  pore  fluid 
flow  in  the  laminar  regime  and  Equation  2-36  describing  flow  in  the  turbulent 
regime.  As  shown  in  Section  6,  our  experimental  work  demonstrated  that  pore 
fluid  transitions  rapidly  into  the  turbulent  regime  under  only  modest  pore 
pressure  gradients.  For  practical  purposes,  in  large  amplitude  dynamic 
loadings,  any  influence  of  the  mass  increment  factor,  r,  is  completely 
obscured  by  the  transition  from  laminar  to  turbulent  flow.  We,  therefore,  can 
use  Equation  2-36  to  represent  flow  in  both  the  laminar  and  turbulent  regimes. 

Substituting  the  flow  parameters  for  Equation  2-36  into  Equation  2-34 

gives 


JM  °  £  *i  ♦  gj  *i2  ♦  ♦  PfUi  (3-4) 

Equation  3-4  can  be  expressed  in  the  following  simple  form  by  using  Equations 
2-19  and  2-22 


(3-5) 


and 


_ k 

k'  =  yf  |  l  +  k3*  |  Wi|  J  (3-6) 

where  k'yf  is  an  equivalent  permeability  applicable  in  the  turbulent  flow 
regime. 

To  simplify  mathematical  derivations,  the  k'  in  Equation  3-6  is  evaluated 
at  time  step  n-1  and  this  value  is  used  at  time  step  n.  k'  is  then  recomputed 
at  time  step  n  for  use  during  the  next  time  step,  n+1. 

For  time  step  n 


k'=  k'n_i 


(3-7) 


At  this  point  we  are  ready  to  express  the  absolute  fluid  velocity  in 
terras  of  the  pore  pressure  gradient  and  skeleton  displacement.  The  apparent 
fluid  velocity,  w-j,  can  be  expressed  in  terras  of  absolute  fluid  velocity,  U^, 
and  skeleton  velocity,  U{,  by  differentiating  Equation  2-21  with  respect  to 
time. 


»  n(U,  -  G,) 


Substitution  of  Equation  3-8  into  Equation  3-5  yields 

PfUi  *  £7  Gi  23  w.  i  * 


(3-8) 


(3-9) 


By  defining  the  above  fluid  flow  equation  at  time  step  n-l$, 

pfUin“*S  +  Jr  <iin_*  =  6<n_*  (3-10) 

The  field  variable  at  time  step  n-Jj  can  be  approximated  by  the  field  variables 
at  time  steps  n-1  and  n  as;  follows; 
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.  L  (6/  -  Ui"-!) 

(3-lla) 

0in-*  =  H  (U^n  +  tj.n-1) 

( 3—1 lb ) 

+  Tr^n-l) 

(3-llc) 

^i0-3*  3  ^  (Au.^) 

(3-lld) 

Substituting  Equation  3-11  into  Equation  3-10  and  solving 
fluid  velocity  at  time  step  n,  we  obtain 

for  the  absolute 

<iin  3  ai  U^-l  +  a2  +  n,^n)  +  a3  Au^n 

(3-12) 

where 


a,  «  2k>pf  "  n  At 
2k 'pf  +  n  At 


a2 


k'  At 


2k'pf  ♦  n  At 


(3-13) 


a3 


2n 


2k'pf  ♦  n  At  f 


Now  combining  Equations  3-llb,  3-llc,  3-lld,  3-12,  and  3-10  and  solving  for 
the  absolute  fluid  acceleration  at  time  step  n-ls,  we  obtain 


PfU^n~*<  =  a4  U^n"l  +  a5  (st,  ♦  rr,^n)  ♦  a6  Au^n 
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where 


a4  =  "2k7  (1+al} 


a6  =  *  (1  -  £r  a2) 

n _  ,,  At  _  . 

ag  =  k'At  (1  “  2  3) 


(3-15) 


At  this  point,  we  have  defined  the  motion  of  the  pore  fluid  in  terms  of 
the  pore  pressure  gradients  and  skeleton  motion  and  can  determine  all  aspects 
of  the  skeleton  and  pore  fluid  behavior.  We  are  now  ready  to  derive  the 
global  equilibrium  equations  for  the  bulk  mixture  and  the  pore  fluid. 

3.5  GLOBAL  EQUILIBRIUM  EQUATION  FOR  THE  .BULK  MEDIUM 

Shown  schematically  in  Figure  3.1  are  the  total  stresses  and  virtual 
displacements  on  the  boundary  of  an  infinitesimal  element.  The  total  stresses 
are  in  equilibrium  with  the  applied  boundary  tractions.  Taking  the  solid 
skeleton  movement  as  the  virtual  displacement,  6ux,  the  internal  and  external 
virtual  work  must  be  equal.  The  internal  virtual  work  at  time  step  n  (t  =  tn) 
is  given  by 


6Wj  »  jjSeWon)  dv 


(3-16) 


Substituting  Equations  3-1  and  3-3  into  Equation  3-16, 


fiWj  =  {fiu}T  (l  [B]T[DeP][B]dv){Au}n 

'  Jv  ' 


(3-17) 


+  (I  j  [8]T(  { 1 }  -  ^  [OeP]  (l} )  <G>dv){A7f}n 

+  l  f  dv 

Jv 

I 

The  external  virtual  work  at  time  step  n  is  given  by 

5WE  =  -J  MTjFjJn  +  |  {6u}T{T}n  ds  (3-18) 

=  j  [N]T(p  -  npf)[N]  dv)){(i}n 

-  (I  |  CN]T  npf{U}n  dvj  +  l  |  [N]T{T}n  ds 

The  absolute  fluid  acceleration  term  in  Equation  3-18  can  be  explicitly 
expressed  using  the  approximate  relationships  in  Equation  3-14. 


I  jyMT  n  Pf  Win  dv  “  iRf lln-lJi  (3-19) 

where 


lRfl}n-lJ*  0  I  {vlN)Tn  (a4  l°}n-2  *  a5  Uff*i}n-2  +  iw»i )n-l )  +  a6  lAu)n-l  j  dv 


(3-20) 


Substituting  Equation  3-19  into  3-18, 
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(3-21) 


3WE  *  {60}T  ~{l  JjN]T  (p-npf)[N]dv)  {a}n 

"  tRf lJn-1%  +  (l  |s [N3T{ T}n  ds) 


Since  internal  and  external  virtual  work  are  equal, 

SWX  =  6We  (3-22) 

Now,  global  equilibrium  equations  for  the  bulk  mixture  are  obtained  by  substi¬ 
tuting  Equations  3-17  and  3-21  into  Equation  3-22. 


[Mt]{d}n  +  [Kt]  {AQ}n  +  [c]  { Aft} n  =  {Pu}n 


(3-23) 


where 

{puln  =  lpln  -  X  jvlBlTK-l}  dv  -  iRfl}n-lli 

[«t]  °  X  JvMT(P  *nPf)tNl  dv 

[KT]  =  l  [  [B]T[0®Pj[B]  dv  (3-24) 

Jv 

CC]  =  I  |  [B]T  |{i}  -  ^  [D®P] |l) |<G>  dv 

(F>n  =  I  f  CN]T{T}n  ds 
is 

3.6  GLOBAL  EQUILIBRIUM  EQUATIONS  FOR  PORE  FLUID 

Figure  3.2  shows  schematically  the  complementary  virtual  pore  pressures 
and  skeleton  and  fluid  velocities  on  the  boundary  of  the  infinitesimal 


59 


element.  The  internal  fluid  movements  relative  to  the  solid  skeleton  are 
compatible  with  the  specified  boundary  flux.  Thus,  taking  the  complementary 
virtual  stresses  as  the  pore  pressure  field  at  a  certain  time,  t,  the  internal 
and  external  complementary  virtual  work  which  is  done  between  time  tn_i  and  tn 
must  be  the  same.  That  is, 


6UZ*  =  SWE*  (3-25) 


where  the  internal  complementary  work  is  given  by 


6Wj*  =  f  [  6jt  n  rr(ep  -  ev)  dv  dt  +  [  (  { 6tt,  , } t{ w)  dv  dt  (3-26) 

Jtn-!  Jv  Jtn_!  Jv 


and  the  external  complementary  work  by 

dWE*  =  f  6n  Q  ds  dt  (3-27) 

J tn-i  h 


Substitution  of  Equations  2-16  and  3-5  into  Equation  3-26  yields 


6Wj*  a  J  n  |  (a  -  ~^~{ljT[DeP]{lj|n  dv 


(3-28) 


fiirT  /{ i}T  -  ^  {ljT[OeP)  {ejj  dv 


♦ 


{6mr. ^ }T  k'{rr,i)  dv 
v 


-  [  pfk'{U) 

Jv 


dv  dt 


The  absolute  fluid  acceleration  terra  in  Equation  3*28  can  be  explicitly 
expressed  using  Equation  3-14  as  follows? 
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f  f  {  ilT  Pfk'  {ii}  dv  dt  a  f  {6^, ^}Tk'pf{u}n-iJj  dv  *  At 
Jtn-l  1  Jv 


(3-29) 


And  the  third  term  in  Equation  3-28  can  be  integrated  with  respect  to  time  as 
below; 


{jr,i}  dv  •  dt  a  f  {Sn^JTk1  |%{Ajr}n  +  {*}„_! 


dv  •  At 
(3-30) 


Now,  substituting  Equations  3-29  and  3-30  into  3-28  and  discretizing  by 
Equation  3-1,  we  have  the  following  expression; 


$*1*  =  [ [e]  {AJt}n  -  [c]T{A0}n 

+  [H]  (*At  {AR}n  ♦  At  (3-31) 

•  I  fvWTk'  pf  (^)n-Ui  * -dv  •  At 


where 

c  ^ 

[E]  a  l  |  <0>T  | a  -  -J-  {l}T(OeP]{lJ  )  <0>  dv  (3-32) 

CM]  a  l  f  [A]T  k'[A]  dv 

Jv 

for  the  external  coopleaentary  virtual  work,  discretizing  Equation  3-27  by 
Equation  3-1  yields 


where 


6We*  *  (6S)T 


(3-33) 


(Q|  =  X  J  <G>T  Q  ds 


(3-34) 
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Assuming  a  linear  variation  of  {Q}  between  time  tn_!  and  tn>  Equation  3-33 
can  be  expressed  as 


SWE*  =  {57f}T 


(3-35) 


Now,  substituting  Equations  3-31  and'3-35  into  Equation  3-25,  the 
following  global  equilibrium  equations  for  the  pore  fluid  are  obtained: 


[C]T{AQ}n  +  |-[e]  -§*[«])  [M)n  =  (P„}n 


(3-36) 


where 

lpn)n  s  At  M’Wn-l  ~  i Rf 2 “  -p  (  ^n-1  +  l^Jn) 
and 

{Rf2Jn-l*»  3  I  |  [A] Tk ' At  •  |  a4ju}n.2 

♦  a5  ((K,i Jn-2  *  lff»i ) n— 1 )  +  a6  IMn-1  )  dv 


(3-37) 


(3-38) 


3.7  COX8INED  GLOBAL  EQUILIBRIUM  EQUATIONS 

Equations  3-23  and  3-36  can  be  combined  in  the  following  matrix  form: 


[H]{3}n  *  [K] {Ad|n  »  {P}n 


(3-39) 


where  the  mass  matrix  [m]  is 


62 


[Mt]  I  [0] 

M  -  - I 


(3-40) 


the  stiffness  oatrix  [k]  is 

1  * 

[*<t] 

1  [c] 

w »  -  -  -  - 

1 . 

(c]T 

the  force  vector  (P}n  is 

» 

lpuln 

|Pln  •  .  _ 

IMn 

the  acceleration  vector  {d)n  is 

(uln 


and  the  displacement  increment  vector  {4d|n  is 

UuJn  I 


Mn 


l^ln 
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3.8  LINEARIZED  GLOBAL  EQUILIBRIUM  EQUATIONS 


Introducing  a  time  integration  method  which  incorporates  both  Newmark's 
8  method  and  Wilson's  6  method,  the  generalized  acceleration  vector  is 
expressed  as 


{d}n  =  CjlAdJn  ♦  C2(d}n-i  *  C3{d} n-l 


where 

1 

1  =  ee3& t* 

C3  a  1  v  200 

and  the  generalized  velocity  vector  is  expressed  as 

|3}n  3  BI  {Ad}n  t  8j  {d}n-i  t  83  {d}n-j 


(3-41) 


(3-42) 


(3-43) 


where 


00  3&t 


Bj  *  1 


.  JL 
002 


(3-44) 


®3  3  1  "  200  6t 

Substituting  Equation  3-41  into  Equation  3-39  and  rearranging,  we  can 
obtain  the  following  linearized  global  equilibrium  equations  which  can  be 
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solved  simultaneously  at  each  step: 


[K]{Ad}n  =  (P}n 


where  the  generalized  stiffness  matrix  [k]  is 

[K]  =  CiCM]  +  [K] 


(3-45) 


(3-46) 


and  the  generalized  force  vector  {P}n  is 

(P}n  =  (Pin  ~  [M]  Cgjd}^!  +  C3{H}n~1 

3.9  SUMMARY  OF  GLOBAL  EQUILIBRIUM  EQUATIONS 

A  summary  of  the  global  equilibrium  equations  includes: 
1.  Bulk  Mixture  Equilibrium 


[Mt]{a}n  +  [Kt]{A0}n  *  [c]{AJt}n  -  {Pu}n 


(3-23) 


2.  Pore  Fluid  Equilibrium 


[C]T|AQ|n  .  -[6]  -  f  [H]  |A»|n  =  (P„)n 


(3-36) 


3.  Combined  Equilibrium  Equation 


[i<]{Ad}n  -  {P}n 


(3-46) 
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Figure  3.1.  Schematic  depiction  of  total  stress  and  virtual  displacement 
field  in  x-direction. 
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SECTION  4 


MPOAP  VERIFICATION  PROBLEMS 


4.1.  INTRODUCTION 

The  formulations  described  in  Section  3  have  been  implemented  in  the  two- 
dimensional  multiphase  finite  element  code  MPDAP.  Features  of  MPOAP  include: 

1.  Special  configurations  to  solve  plane  strain,  axi-symmetric  and 
spherically  symmetric  problems; 

2.  Both  decoupled  and  fully  coupled  compressibility  models.  The  simpler 
decoupled  model  treats  the  mixture  and  skeleton  compressibilities 
independently  (analogous  to  parallel  springs).  The  fully  coupled 
model  includes  the  influence  of  grain  compression  due  to  pore  pressure 
and  the  influence  of  effective  stresses  on  the  compressibility  of  the 
skeleton  (see  Blouin  and  Kim,  1984). 

3.  Nonlinear  pore  water  and  solid  grain  compressibility  models  (Kio, 
Blouin  and  Timian,  KC.). 

4.  A  variety  of  available  material  models  describing  the  drained  skeleton 
response.  These  include: 

•  linear  elastic; 

•  decoupled  nonlinear  material  model  which  treats  volumetric  and 
deviatoric  shear  response  as  independent  (i.e.  no  di latency) 

(  Kim,  Blouin  and  Timian,  1966); 

•  elasto-plastic  material  model  with  a  nonlinear  strength  envelope 
(models  di latency  with  an  associated  flow  rule  (Kio,  Blouin  and 
Timian,  1986);  and 

•  ARA  three  invarient  plasticity  model  with  strain  hardening  and 
softening  (models  di latency  with  a  nonassociated  flow  rule) 

(Herkle  and  Oass,  1986). 

5.  Single  phase,  two  phase  (fully  saturated)  and  three  phase  (partially 
saturated)  capability; 
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6.  Static,  dynamic  and  quasi-static  (consolidation)  capability; 

7.  Calculation  efficiency  features  including; 

•  a  single  point  stress  calculational  technique  for  rapid  computation 
of  stresses,  strains,  and  nonlinear  constitutive  properties  (Kira, 
Blouin  and  Timian,  1986);  and 

•  use  of  pore  pressure  at  nodal  points  to  efficiently  compute  relative 
pore  fluid  motions; 

8.  Nonlinear  fluid  friction  model  for  both  laminar  and  turbulent  pore 
fluid  flow  conditions. 

A  users  manual  for  MPDAP  is  included  in  Appendix  A,  This  describes 
available  options  and  input  format. 

In  this  section,  we  describe  four  problems  used  in  the  initial  verifica¬ 
tion  of  MPOAP  to  check  both  the  theoretical  formulations  and  their  numerical 
implementation.  These  include: 

1.  an  undrained  compressibility  problem; 

2.  spherical  wave  propagation  in  an  elastic  single  phase  material; 

3,  a  linear  quasi-static  consolidation  problem; 

4,  two  phase  wave  propagation  in  saturated  soil  and  rock  having  linear 
elastic  skeletons. 

In  these  problems,  the  HPDAP  calculations  are  compared  with  closed  form 
solutions  and/or  numerical  solutions  generated  with  the  code  TPOAPII.  TPOAPII 
is  our  previous  two  phase  code,  described  by  Kim,  Blouin  and  Timian  (1986), 
which  does  not  have  most  of  the  advanced  features  incorporated  into  HPOAP. 
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4.2.  VERIFICATION  PROBLEM  1,  CONSTRAINED  COMPRESSIBILITY  OF  UNORAINED 

SATURATED  GRANULAR  SOILS. 

The  first  verification  problem  is  to  demonstrate  analytically  that  the 
MPOAP  formulations  in  Section  3  degenerate  to  the  fully  coupled  solution 
(Blouin  and  Kim,  1984)  for  undrained  uniaxial  strain  loadings  as  schematically 
shown  in  Figure  4.1. 

To  simplify  the  analytical  derivation,  three  degrees  of  freedom  are 
introduced;  the  first  degree  of  freedom  represents  the  solid  skeleton 
displacement  and  the  second  and  third  degrees  of  freedom  represent  pore  fluid 
pressures.  Under  uniaxial  strain  static  undrained  loading  conditions,  the 
general  MPDAP  formulation  (Equation  3-39)  degenerates  to  the  following  simple 
equation: 


(4-1) 


For  the  single  element  with  unit  length.  Equation  4-1  can  be  expressed  as 
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(4-4) 
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1_  .  Ks 
Km  Kg^ 


Ms  =  constrained  modulus  of  skeleton 
Uj  =  skeleton  displacement 
jt o  =  pore  fluid  pressure 
ov  =  applied  total  stress 

Solving  Equation  4-2  for  the  pore  fluid  pressure,  we  have 


^0  -  ov 


M 


S 


Substitution  of  Equation  4-3  and  4-4  into  4-5  yields 
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where 


Kg2Ms  ♦  KmKs2 
-o  _ _ 
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(4-7) 


Equation  4-6  is  identical  to  the  closed  form  solution  derived  by  Blouin  and 
Kim  (1984). 

4.3.  VERIFICATION  PROBLEM  2,  ELASTIC  SPHERICAL  WAVE  PROPAGATION  IN  A  SINGLE 
PHASE  MEDIUM. 

The  purpose  of  this  verification  problem  is  to  check  the  global 
equilibrium  equations  for  the  bulk  mixture  (Equation  3-23)  in  the  one- 
dimensional  spherical  coordinate  system  when  the  pore  fluid  is  not  present. 

Figure  4.2  shows  a  12  inch  hollow  spherical  hole  in  an  infinite  elastic 
medium  subjected  to  a  100  psi  internal  <tep  load.  Material  properties  and 
tine-steps  used  for  the  calculations  are  included  in  Figure  4.2.  In  Figure 
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4.3,  radial  stress  profiles  at  5.5  msec  are  plotted.  The  MPDAP  calculations 
give  good  agreement  with  the  closed  form  solution  except  for  some  smearing  of 
stresses  at  the  shock  wave  front.  This  smearing  could  be  minimized  by 
reducing  the  element  and  time  step  sizes. 

4.4.  VERIFICATION  PROBLEM  3,  ONE-DIMENSIONAL  LINEAR  CONSOLIDATION. 

This  verification  problem  is  used  to  check  the  quasi-static  flow  portion 
of  the  MPDAP  formulations  against  Terzaghi's  closed  form  consolidation  solu¬ 
tion. 


A  fully  saturated  soil  deposit  is  assumed  to  overlay  a  rigid  impermeable 
base.  The  soil  deposit  is  subjected  to  a  step  loading  at  the  ground  surface, 
which  is  assumed  to  be  a  free-draining  boundary.  The  initial  excess  pore 
pressure  distribution  is  assumed  to  be  constant  throughout  the  deposit. 

Twenty  equally  spaced  4-node  elements  are  used  with  a  non-dimensional  time 
increment  factor  AT  a  0.005.  Plotted  in  Figure  4.4  is  the  profile  of  the  nor¬ 
malized  excess  pore  water  pressure  at  time  factor  T  =>  0.5,  at  which  time 
about  75$  of  the  excess  pore  pressure  is  dissipated  on  average  throughout  the 
soil  deposit.  The  calculated  excess  pore  water  pressures  show  close  agreement 
with  Terzaghi's  exact  solution. 

4.5.  VERIFICATION  PROBLEM  4,  ONE-DIMENSIONAL  WAVE  PROPAGATION  IN  SATURATED 

LINEAR  ELASTIC  SOILS  AND  ROCKS. 

Two  series  of  one-dimensional  calculations  of  a  vertically  propagating 
planar  compression  wave  were  performed  using  MPOAP  on  idealized  saturated  soil 
and  rock.  The  input  loading,  as  shown  in  Figure  4.5,  was  a  short  rise  time 
triangular  pulse  with  a  peak  stress  of  5,000  psi  and  a  positive  phase  duration 
of  10  msec.  The  loading  pulse  was  applied  to  the  saturated  sand  and  rock 
having  the  properties  listed  in  Figure  4.5.  The  load  was  applied  to  an  imper¬ 
meable  boundary  at  the  ground  surface. 

Three  permeability  values  were  used  in  each  series  of  calculations  on  the 
two  materials,  0,001,  0.1,  and  1.0  in/s.  Results  of  the  MPOAP  calculations 
are  compared  to  identical  calculations  using  TPOAPII  in  Figures  4.6  through 
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4.8  for  the  rock  and  4.9  through  4.11  for  the  soil.  Each  series  of  calcula¬ 
tions  compares  pore  pressure  and  effective  stress  profiles  from  the  TPDAPII 
and  MPOAP  calculations  at  times  of  10  msec  and  20  msec  after  application  of 
the  boundary  stress.  Each  figure  {parts  a  through  d)  shows  the  profiles  for 
one  of  the  three  assumed  permeabilities.  The  decoupled  compressibility  option 
in  MPDAP  was  used  because  TPDAPII  does  not  have  a  fully  coupled  compressibi¬ 
lity  model. 

Overall,  the  comparisons  between  the  two  sets  of  calculations  is 
excellent.  As  would  be  expected,  effective  stresses  in  the  soil  are  very 
small  compared  to  the  pore  pressures  because  of  the  high  compressibility  of 
the  soil  skeleton.  The  effective  stresses  in  the  rock  are  higher  than  the 
pore  pressures  because  the  rock  skeleton  is  stiffer  than  the  solid  grain-pore 
water  mixture. 

The  only  major  discrepancy  occurred  in  the  highest  permeability  soil 
calculations  shown  in  Figure  4.11.  In  the  HPOAP  calculation,  both  pore 
pressure  and  effective  stress  profiles  showed  high  frequency  oscillations. 
These  are  believed  to  occur  because  of  the  high  contrast  between  the  soft  ske¬ 
leton  stiffness  and  the  stiff  bulk  mixture  modulus.  The  high  permeability,  in 
effect,  creates  an  even  higher  contrast  between  these  two.  Further  analysis 
will  be  required  to  better  define  and  correct  this  problem.  For  practical 
values  of  permeabilities  {less  than  0.1  in/s),  the  oscillations  are  not  a 
problem. 

The  high  frequency  oscillations  in  some  of  the  calculations  at  the  imper¬ 
meable  ground  surface  may  be  associated  with  waves  of  the  second  kind.  They 
occur  in  both  the  HPOAP  and  TPDAPII  calculations  and  are  most  noticeable  at 
the  highest  permeabi 1 ity,  as  would  be  expected  for  waves  of  the  second  kind. 


.  A  single  element  subjected  to  undrained  uniaxial  strain  loading. 
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Figure  4.3.  Radial  stress  profile  at  5.5  msec,  verification  problem  7 
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Figure  4.5.  Loading  time  history  and  material  properties  used  in  1-D  plane  strain  elastic  wave 
propagation  through  saturated  soil  and  rock. 
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Figure  4.5a.  Pore  Fluid  pressure  profile  at  10  msec,  rock  with  k  =  0.001  in/sec. 
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Effective  vertical  stress 
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Effective  stress  profile  at  20  msec,  rock  with  k  =  0.001  in/sec 
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Figure  4.7a.  Pore  fluid  pressure  profile  at  10  msec,  rock  with  k  -  0.1  in/sec. 
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Figure  4.7c.  Pore  fluid  pressure  profile  at  20  msec,  rock  with  k 
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Figure  4.7d.  Effective  stress  profile  at  20  msec,  rock  with  k  =  0.1  in/sec. 
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Figure  4.8a.  Pore  fluid  pressure  profile  at  10  msec,  rock  with  k  =  1.0  in/sec. 
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Figure  4.8b.  Effective  stress  profile  at  10  msec,  roc 
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Effective  stress  profile  at  20  msec,  rock  with  k 
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Figure  4.9b.  Effective  stress  profile 
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4.9c.  Pore  fluid  pressure  profile  at  20  msec,  soil  with  k  =  0.001  in/sec 
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Figure  4.9d.  Effective  stress  profile  at  20  msec,  soil  with  k  =  0.001  in/sec. 
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pressure  profile  at  20  msec,  soil  with  k  =  0.1  in/sec. 
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TWO-PHASE  DYNAMIC  ANALYSIS  (  MPDAP  VS 
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Figure  4.11b.  Effective  stress  profile  at  10  msec,  soil  with  k  =  1.0  in/sec. 
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Effective  stress  profile  at  20  msec*  soil  with  k  K  1,0  in/sec. 
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SECTION  5 


NUMERICAL  AND  THEORETICAL  TREATMENT  OF  WAVES 
OF  THE  FIRST  AND  SECOND  KIND 


5.1.  INTRODUCTION 

Kim,  Blouin  and  Timian  (1987)  presented  both  a  theoretical  and  numerical 
treatment  of  wave  propagation  and  damping  in  saturated  porous  media.  In 
Section  4  of  that  report,  a  closed  form  solution  for  wave  propagation  velocity 
and  damping  in  fully  saturated  porous  media  having  elastic  skeleton  properties 
was  derived  for  a  fully  coupled  model  with  linearly  compressible  solid  grains 
and  pore  water.  This  solution  demonstrated  the  existence  of  two  types  of 
compression  waves,  termed  waves  of  the  first  and  second  kinds.  It  also  pro¬ 
vided  a  means  of  benchmarking  and  verifying  multiphase  code  calculations  and  a 
means  of  planning  and  guiding  a  further  investigation  of  these  two  types  of 
waves  using  our  finite  element  two-phase  codes. 

During  the  past  year,  we  have  expanded  on  this  initial  study.  We  per¬ 
formed  a  numerical  investigation  of  waves  of  the  second  kind  in  saturated 
porous  rock  to  determine  some  of  the  characteristic  properties  of  waves  of 
the  second  kind.  In  addition,  the  method  of  characteristics  was  used  to 
verify  and  expand  our  previous  theoretical  derivation  for  a  general,  rather 
than  a  harmonic,  loading  function.  Results  of  these  numerical  and  theoretical 
studies  are  described  in  this  section. 

5.2.  NUMERICAL  CALCULATIONS  OF  WAVES  OF  THE  FIRST  AND  SECOND  KINO 

The  general  theoretical  solution  for  compressional  wavespeeds  and  damping 
was  incorporated  into  the  numerical  code  TWAVE,  described  by  Kim,  Blouin  and 
Timian  (1987).  TWAVE  was  used  in  a  parametric  study  of  the  influence  of  exci¬ 
tation  frequency  and  variations  in  material  properties  on  propagation  velocity 
and  damping.  Compressional  wave  velocity  for  waves  of  the  first  kind,  as 
shown  in  Figure  5.1,  was  found  to  vary  as  a  function  of  the  frequency- 
permeability  product,  with  a  zone  where  wavespeed  transitions  from  a  lower 
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bound  value  to  a  higher  bound  value  with  increasing  values  of  the  product. 
Damping  is  seen  to  be  at  a  maximum  where  the  rate  of  change  in  wavespeed  is 
greatest.  As  shown  in  Figure  5.2,  waves  of  the  second  kind  also  have  a  tran¬ 
sition  in  wavespeed  from  near  zero  at  low  values  of  the  frequency-permeability 
product  to  an  upper  bound  value  at  higher  values  of  the  product. 

These  theoretical  predictions  were  also  observed  in  a  set  of  numerical 
calculations  performed  by  TPOAPII.  A  series  of  3  one-dimensional  calcula¬ 
tions,  similar  to  those  described  in  Section  4.5,  for  a  vertically  propagation 
planar  compression  wave  were  performed  for  three  different  permeabilities. 
Permeabilities  were  selected  sc  that  the  frequency-permeability  product  fell 
in  the  lower  bound,  upper  bound  and  transitional  wavespeed  regimes.  The  input 
loading  was  the  same  as  described  in  Figure  4.5  and  consisted  of  a  short  rise 
time  triangular  pulsa  with  a  peak  stress  of  5,000  psi  and  a  positive  phase 
duration  of  10  msec.  The  loading  pulse  was  applied  to  saturated  soil  having 
the  properties  also  listed  in  Figure  4.5.  Snapshots  of  the  pore  water 
pressure  profiles  at  four  different  times,  from  10  to  40  msec,  rre  shown  in 
Figure  5.3,  Three  calculations  are  shown  for  permeabilities  of  0.001,  0.1  and 
10  in/s. 

The  trends  in  these  calculations  are  similar  to  those  predicted  from  the 
TWAVE  closed  form  solution.  The  wavespeed  increases  substantially  with 
increasing  permeability.  For  the  lowest  permeability,  the  velocity  of  the 
wavefront  is  about  5300  ft/s  and  for  the  highest  permeability,  it  is  about 
5900  ft/s.  The  wavespeed  computed  for  this  material  according  to  the  decoupled 
undrained  modulus  described  by  Blouin  and  Kim  (1984)  is  about  5200  ft/s  which 
agrees  very  well  with  the  low  permeability  calculation. 

There  is  a  dramatic  alteration  of  the  wave  shape  for  the  intermediate 
permeability  calculation.  The  wavefront  is  smeared  and  the  amplitude  is 
significantly  attenuated  relative  to  the  calculations  in  the  less  permeable 
and  in  the  more  permeable  materials.  This  is  a  clear  indication  that  excess 
damping  occurs  in  the  transition  region  between  the  lower  bound  and  upper 
bound  wavespeeds,  as  predicted  by  the  TWAVE  closed  form  solution. 
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An  additional  calculation  was  performed  using  TPDAPII  which  was  tailored 
to  highlight  the  phenomenology  associated  with  compressional  waves  of  the 
second  kind.  In  order  to  accentuate  the  wave  of  the  second  kind,  a  permeabie 
rock  skeleton  was  used  having  properties  listed  in  Figure  5.4.  The  saturated 
rock  was  loaded  with  a  downward  sinusoidal  pressure  pulse  having  a  peak 
pressure  of  5000  psi  and  a  frequency  of  50  Hz  and  a  half  sine  duration. 

Figure  5.4a  (Kim  and  Blouin)  shows  the  pore  pressure  profiles  for  both  the 
waves  of  the  first  and  second  kinds  at  20  msec.  The  pore  water  is  in 
compression  in  both  instances,  with  the  lower  amplitude  wave  of  the  second 
kind  lagging  well  behind  the  wave  of  the  first  kind.  The  effective  stress  in 
the  rock  skeleton,  shown  in  Figure  5.4b  is  in  compression  in  the  wave  of  the 
first  kind,  but  in  tension  in  the  wave  of  the  second  kind.  The  combined  wave¬ 
forms  are  shown  in  Figure  5.4e  where  it  is  obvious  that  the  compressional 
stresses  in  the  pore  water  of  the  second  wave  are  balanced  by  the  tensile 
stresses  in  the  porous  skeleton. 

Examination  of  the  TPDAPII  output  for  the  pore  water  and  skeleton  motion 
shows  that  in  the  wave  of  the  first  kind,  both  the  pore  water  and  skeleton  are 
moving  in  the  direction  of  the  wave  propagation  and  are  both  in  compression. 
However,  in  the  wave  of  the  second  kind,  the  motion  of  the  skeleton  and  pore 
water  are  out  of  phase.  The  pore  water  is  in  compression  and  is  moving  in  the 
direction  of  the  wave  propagation.  The  skeleton  is  in  tension  and  is  moving 
opposite  of  the  direction  of  propagation.  Thus,  waves  of  the  second  kind 
appear  to  consist  of  a  surge  of  pore  water  moving  through  the  porous  skeleton. 
It  would  be  very  desirable  to  obtain  experimental  verification  of  the  waves  of 
the  second  kind,  as  well  as  verification  of  the  dependence  of  wavespeed  for 
waves  of  the  first  kind  on  the  frequency-permeability  product. 

5.3.  THEORETICAL  DERIVATION  Of  FIRST  AND  SECOND  KIND  KAVESPEEDS  6V  THE  METHOD 

OF  CHARACTERIZATION 

5.3.1.  Introduction 

Kim,  Blouin  and  Tiraian's  (1997)  derivation  of  wavespeed  and  damping  for 
waves  of  the  first  and  second  kind  utilized  a  harmonic  excitation  function. 

In  order  to  further  verify  the  theoretical  existence  of  waves  of  the  second 
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kind,  we  have  developed  an  alternate  derivation  using  the  method  of  charac¬ 
teristics  which  doesn't  rely  on  a  harmonic  excitation  function.  This  solution 
confirms  the  existence  of  waves  of  the  second  kind  in  saturated  porous  media 
with  no  fluid  friction  damping.  The  wavespeeds  obtained  in  this  solution  are 
identical  to  the  upper  bound  wavespeeds  from  the  more  general  1987  solution 
where  the  frequency-permeability  product  tends  toward  infinity. 

5.3.2.  Derivation  of  Wave  Propagation  in  Saturated  Media 

The  four  governing  field  equations  for  uniaxial  strain  loading  from  Kim, 
Blouin  and  Timian  (1987)  (Equations  4-35,  4-36,  4-40  and  4-46)  are: 


whe  governing  differential  equilibrium  equation  for  the  bulk  mixture: 


6W; 


Pf 


at 


a<x2 

32 


=  0 


(5-1) 


the  governing  differential  equation  for  motion  of  the  pore  fluid: 


8uz  if  Pf  3l!iz  3rr 

“for  1  r  *2  *  ~  »*r|  sr  -  « * 0 


the  governing  differential  equation  for  continuity  of  flow: 


U-CgKs) 


(5-3) 


the  governing  differential  equation  for  the  effective  stress  law: 


3u2  arr  ?<rz 

"ssr  *  <l-Vs>at  -  aT  “  0  (S“" 

NBi  a  and  n  are  both  positive  in  tension 

Equation  5-1  through  5-4  are  four  linear,  first  order,  partial  differential 
equations  involving  the  quantities  u2.  wz,  oz,  and  n,  each  of  which  has  par- 
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tial  derivatives  with  respect  to  z  and  t.  Since  we  have  eight  partial  deriva¬ 
tives,  but  only  four  equations,  we  need  four  more  equations,  which  are: 


3uz  3uz 

*  aT  d2  *  5T  dt 


(5-5) 


3wz  3wz 

d*z  ’  5T  dz  *  aT  dt 


(5-6) 


doz 


dz  ♦ 


(5-7) 


and 


art 

dn  s  az  dz  * 


(5-8) 


Equations  5-1  through  5-8  can  be  written  in  matrix  form  as  follows: 


(5-9) 


Equation  5-9  provides  a  means  of  determining  the  eight  partial  derivatives, 


^  /  W2  2  ^  ^  2  /  ^2  ^  >  7T  2 

at  the  point  and  the  differentials  du2, 


and  7T  t  at  a  point,  when  w2  is  known 
dw2,  da2  and  dn  are  specified  in  a 


given  direction 


dz  =  cdt  (5-10) 

except  when  the  value  of  c  renders  the  coefficient  matrix  singular. 

For  convenience,  let 


£l(l+r)  =  a 


(5-11) 


1  -  CgKs  =  b 


(5-12) 


Cq2Ks  -  ~  =  d 


(5-13) 


then  the  coefficient  matrix  in  Equation  6-9  is 


0 

0 

b 

Hs 

cdt 

0 

0 

0 


P 

Pf 

0 

0 

dt 

0 

0 

0 


0 

0 

1 

0 

0 

cdt 

0 

0 


Pf 

a 

0 

0 

0 

dt 

0 

0 


-1 

0 

0 

0 

0 

0 

cdt 

0 


0 

0 

0 

-1 

0 

0 

dt 

0 


0 

-1 

0 

0 

0 

0 

0 

cdt 


0 

0 

d 

b 

0 

0 

0 

dt 


(5-14) 
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Multiplying  each  even  column  of  |  A  J  by  c,  and  subtracting  it  from  the  previous 
column  leaves  the  values  of  I  A  I  unaltered,  so  that 


P  -cpf 
Pf  -ca 


I A I  =  A  = 


dt  0 
0  0 


(5-15) 


Using  the  elements  dt  in  the  5th  through-. 8th  rows  as  the  pivotal  elements  in  a 
Laplace  expansion  yields 


A  a  (dt)4 


!+cpf  +ca 
(+1)  b  1 


-CP  -cpf 

-1 

0 

-cpf  -ca 

0 

-1 

b  1 

0 

-cd 

Ms  0 

c 

-cb 

i  +1 

+ep 

♦Cpf  0 

-cd  +  (c) 

-CPf 

-ca  -1 

-cb 

b 

1  -c< 
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s 


(-c2pfb  -  c2adMs  +  c2ab2  -  Ms) 


<dt)4 


+  c(c3pad  -  cpfb  +  cp  -  c3pf2d) 


=  (dt)4 


(pad  -  Pf2d)c4  +  (-pfb  -  adMs  +  ab2  -  pfb  +  p)  c2  -  Ms 


=  -(dt)4  (Pf2d  -  pad)c4  -  (ab2  +  p  -  2pfb  -  adMs)  c2  +  Ms 


(5-16) 


If  we  set 


pf2d  -  pad  =  -d(pa  -  pf2)  =  A' 


(5-17) 


ab2  +  p  -  2pfb  -  adMs  =  2B 


(5-18) 


Me  =  C 


(5-19) 


then  the  condition  that  A  vanish,  which  renders  A  singular,  reduces  to: 


A'c4  -  2Bc2  +  C  =  0 


(5-20) 


which  means  that 


2  B  t  /B2  -  4A'C 
3  2A' 


Cc  3 


(5-21) 


This  wavespeed  solution  is  identical  to  the  upper  bound  wavespeed  derived  by 
Kim,  Blouin  and  Timian  (1987)  given  in  Equation  4-111.  The  plus  sign  produces 
upper  bound  velocities  of  waves  of  the  first  kind  in  an  undamped  medium  and 
the  minus  sign  produces  upper  bound  velocities  of  waves  of  the  second  kind  in 
an  undamped  medium. 
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As  compatibility  checks,  we  consider  two  cases:  where  there  is  no  pore 
fluid  and  when  there  is  no  skeleton  (water  only).  When  there  is  no  pore  fluid 


pf  =  0  a  =  0 

r  =  0  b  =  1 

Cg  =  0  d  =  0 

P  =  (l-n)pg 

and  Equation  £-2  reduces  to 


pc 


2  _ 


=  0 


so  that 


c 


A'  =  0 
2B  =  ps 
C  =  Ms 


(5-22) 


(5-23) 


which  is  the  compress ional  wavespeed  in  a  single  phase  dry  medium. 
When  there  is  no  skeleton 

n  =  1  A*  a  o 


a  =  pf 

b  =  1 


2B 


P  f  M 
Kf 


c  a  Ms 


and  Equation  5-20  reduces  to 


Ms  =  0 


(5-24) 
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so  that,  besides  the  trivial  solution  Ms  =  0, 


C2 

Pf 


(5-25) 


which  is  the  compress ional  wavespeed  in  a  pure  fluid. 
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Figare  5.3.  Influence  of  permeability  on  wavespeed  and  damping  for  saturated  sand 
with  constant  skeleton  constrained  modulus  of  6000  psi  and  porosity  of 
351  (Kira,  Blouin  and  Trmian,  1987). 


SECTION  6 


EXPERIMENTAL  EVALUATION  OF  FLUID  FLOW  IN  DUCTS  AND  SOILS 

6.1.  INTRODUCTION 

As  described  in  Section  2.2.2,  the  equation  of  motion  of  the  pore  fluid 
(Equation  2.34)  is  one  of  the  four  field  equations  required  to  describe 
multiphase  dynamic  response.  Nearly  all  previous  analysis  in  this  subject 
area  was  based  on,  or  heavily  influenced  by,  Biot's  pioneering  theoretical 
work  on  laminar  flow  (1956,  1962a,  1962b).  The  experimental  findings  pre¬ 
sented  in  this  section  warrant  a  move  away  from  Biot's  formulation  of  fluid 
friction,  simply  because  turbulent  flow  is  found  to  govern  the  fluid  friction 
equation  for  the  range  of  pore  pressure  gradients  generated  by  explosive  or 
rapid  dynamic  loadings.  The  consequences  of  these  experimental  findings  on 
the  theoretical  work  and  its  implementation  in  MPQAP  are  described  in  Sections 
2  and  3. 

Two  sets  of  experimental  results  are  presented  in  this  section:  a  sum¬ 
mary  of  the  steady  state  and  dynamic  test  results  in  idealized  flat  ducts  (the 
experimental  analog  to  Biot's  theoretical  formulation)  and  a  set  of  steady 
state  flow  data  for  uniform  grained  sand  over  a  wide  range  of  flow  velocities 
and  pore  pressure  gradients. 

6.2.  DYNAMIC  FLOW  TESTS  IN  FLAT  DUCTS 

In  order  to  conduct  the  dynamic  flow  tests  on  the  ducts  and  granular 
soils,  a  device  was  fabricated  for  conducting  high  pressure  steady  state  and 
dynamic  flow  tests  on  both  ducts  and  soils.  This  device,  shown  schematically 
in  Figure  6.1,  is  capable  of  forcing  pore  fluid  through  a  soil  sample  or  duct 
with  a  pressure  difference  of  up  to  5000  psi.  Originally,  the  device  was 
designed  to  be  operated  by  a  servo-controlled  pump  which  could  produce  flow 
through  the  sample  at  a  controlled  rate  or  change  in  rate.  During  the  past 
year,  the  device  was  modified  to  utilize  high  pressure  nitrogen  to  propel  the 
piston  which  forces  the  pore  fluid  through  the  test  specimen.  In  this  way, 
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higher  flow  velocities  could  be  generated  and  very  rapid  changes  in  flow  could 
be  created.  Through  use  of  a  fast  acting  valve,  pressure  gradient  rise 
times  as  short  as  10  msec  could  be  generated.  Additional  details  of  this 
device  are  presented  by  Kim,  Blouin  and  Timian  (1987). 

Using  the  gas  powered  modified  permeability  device,  an  additional  series 
of  11  steady  state  flow  tests  were  run  at  high  flow  velocities  of  up  to  3700 
in/s  in  the  flat  duct.  The  flat  duct  has  a  height  of  0.04  inches,  a  width  of 
0.28  inches  and  a  length  of  8.71  inches.  These  data  are  shown  in  Figure  6.2 
as  pressure  gradient  normalized  by  the  steady  state  flow  velocity  plotted  ver¬ 
sus  the  steady  state  flow  velocity.  They  have  been  combined  with  the  lower 
velocity  data  from  Kim,  Blouin  and  Timian  (1987)  which  were  mostly  run  at  less 
than  1800  in/s  velocity.  The  original  data  fit  is  also  a  good  match  to  the 
higher  velocity  data.  The  intercept  of  the  data  fit  (a  =  .0025)  can  be  used 
to  compute  the  coefficient  of  permeability,  k,  according  to  Equation  2.36  as 

k  =  p.  (6-1) 

a 

as  described  by  Kim,  Blouin  and  Timian  (1987).  A  value  of  k  for  kerosene  used 
in  these  tests  is  11.85  in/s.  This  compares  closely  to  the  analytical  value 
of  11.92  ra/s  computed  in  the  above  report. 

The  slope  of  the  fit  to  the  data  in  Figure  6.2,  b,  can  be  used  to  compute 
the  turbulent  fluid  friction  coefficient,  from  Equation  2-36  as 

/3f  =  b/k“  (6-2) 

For  the  measure  slope  of  1.5  x  10“5  Ib-s^/in®,  a  value  of  =  5.16  x  10"6 
lb-sl,5/in4,5  is  obtained. 

Thus,  the  steady  state  flow  tests  provide  us  with  the  first  two  coef¬ 
ficients  jf  the  general  dynamic  flow  Equation  3-4  which  can  be  rewritten  as 
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Frictional 

Component 


Inertial 

Component 


In  order  to  completely  validate  the  above  flow  equation,  a  series  of  12 
dynamic  flow  tests  were  conducted  on  the  flat  duct.  Representative  data  from 
five  of  these  tests  are  presented  in  Figures  6.3  through  6.7.  During  these 
dynamic  tests,  a  pore  pressure  gradient  is  rapidly  applied  by  the  piston 
driven  by  the  compressed  nitrogen.  Both  the  pore  pressure  gradient  and  flow 
are  monitored  as  functions  of  time  as  shown  in  the  measured  pressure  gradient 
and  flow  displacement  plots  of  part  a  of  Figures  6.3  through  6.7. 

From  the  measured  displacement,  the  fluid  velocity  and  acceleration  are 
obtained  by  differentiating  with  respect  to  time,  as  shown  in  the  bottom  plot 
of  part  b  of  Figures  6.3  through  6.7.  The  velocity  profile  is  then  substi¬ 
tuted  into  Equation  6-3  to  compute  the  frictional  component  of  the  pore 
pressure  gradient  using  the  previously  measured  friction  coefficients,  a  and  b 
of  Equations  6-1  and  6-2.  This  result  is  expressed  in  normalized  form  as  the 
frictional  component  of  the  pressure  gradient  in  the  top  plot  of  part  b  of 
Figures  6.3  through  6.7. 

Using  the  measured  acceleration  from  Figures  6.3  through  6.7,  part  b,  the 
inertia  force  from  the  last  term  of  Equation  6-3  is  computed.  The  normalized 
inertial  component  plus  the  normalized  friction  component  combine  to  equal  the 
normalized  pressure  gradient  of  unity  shown  in  the  plot  of  Figures  6.3  through 
6.7,  part  b.  The  total  normalized  gradients  computed  from  the  velocity  and 
acceleration  data  compare  to  within  about  20%  with  the  actual  measured 
pressure  gradients,  providing  experimental  verification  of  our  dynamic  pore 
fluid  flow  equation  (Equations  6-3  and  3-4). 

Both  the  measured  pressure  gradient  and  computed  frictional  component  of 
the  pressure  gradient  are  normalized  by  the  measured  pressure  gradient  com¬ 
puted  from  the  velocity  and  acceleration  data  in  Figure  6.3  through  6.7.  Note 
that  during  the  high  early  time  acceleration,  the  pressure  gradient  is 


dominated  by  the  inertial  component,  and  as  accelerations  diminish,  the  fric¬ 
tional  component  (dominated  by  the  second  term  of  Equation  6-3)  governs.  The 
agreement  between  the  measured  pressure  gradient  and  that  computed  from  the 
velocity  and  acceleration  data  is  generally  within  20%  throughout  the  entire 
test. 


6.3.  FLOW  TESTS  IN  SATURATED  UNIFORM  GRAINED  SOIL 

A  series  of  flow  tests  was  conducted  on  a  uniform  carbonate  sand  which 
duplicated  the  test  procedures  used  for  the  flat  duct  flow  tests.  The  objec¬ 
tives  of  these  tests  were  to: 

1.  measure  the  fluid  friction  flow  parameters  in  actual  soil  samples; 

2.  determine  where  the  transition  from  laminar  to  turbulent  flow  occurs 
in  a  uniform  sand; 

3.  determine  whether  Biot's  frequency  dependent  laminar  flow  parameters 
are  a  significant  factor  in  determining  the  fluid  friction  in  uniform 
sand;  and 

4.  validate  the  pore  fluid  flow  equation  (Equation  6-3)  in  saturated 
soil . 

In  order  to  make  the  appropriate  fluid  flow  measurements  in  soil,  the 
flat  duct  was  replaced  by  a  cylindrical  soil  specimen  container  having  a 
diameter  of  0.434  inches  and  a  length  of  7.8  inches.  A  uniform  grained  car¬ 
bonate  beach  sand  (Enewetak  beach  sand)  was  obtained  by  sieving  so  that  all  of 
the  pretest  material  was  retained  on  a  number  40  sieve  and  had  a  grain  size 

ranging  from  0.425  to  0.60  mm.  This  is  the  same  sand  used  in  the  grain 

crushing  experiments  described  by  Kim,  Blouin  and  Timian  (1987),  Section  7. 

The  sample  specimens  had  the  following  average  properties: 

dry  density,  =  106.0  lb/ft3 

specific  gravity,  Gs  »  2.81  gm/cm3 

porosity,  n  =  39.5% 
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A  series  of  33  steady  state  flow  tests  were  run  on  the  uniform  beach  sand 
at  apparent  flow  velocities  ranging  fMom  .4  in/s  to  100  in/s  and  pressure  gra¬ 
dients  ranging  from  less  than  1  psi/inch  to  about  250  psi/incn.  Data  from 
these  tests  are  presented  in  Appendix  B  as  apparent  flow  displacement  and 
pressure  gradient  both  plotted  as  functions  of  time.  These  data  are  sum¬ 
marized  in  Table  6.1,  where  the  pressure  gradient  is  also  normalized  by  the 
absolute  flow  velocity.  Figure  6.8  is  a  plot  of  normalized  pressure  gradient 
as  a  function  of  the  absolute  flow  velocity.  The  flow  data  are  conveniently 
fit  by  a  linear  relationship,  indicating  that  flow  is  in  the  turbulent  regime 
over  essentially  the  entire  range  of  pore  pressure  gradients.  Values  of  the 
intercept,  a,  of  Equation  6-1  and  the  slope,  b,  of  Equation  6-2  were: 


a  =  1.380  lb-s/in4 
b  =  .2018  lb-s2/in5 

From  Equation  6-1,  with  yf  for  kerosene  of  0.0291  lb/in3 

k  a  .0211  in/s 


and  from  Equation  6-2, 

jSf  =  .0293  lb«3l‘5/in4*6 

As  a  check  against  the  steady  state  test  data  a  conventional  constant 
head  permeability  test  using  water  as  the  pore  fluid  was  conducted  on  a 
cylindrical  Enewetak  beach  sand  sample  with  a  length  of  6.35  inches  and 
diameter  of  2.50  inches.  A  value  of  k  =  0.0383  in/s  was  obtained  from  this 
test.  In  order  to  convert  the  steady  state  flow  data  to  the  equivalent  per 
meabilility  for  water,  the  following  equation  is  applied. 


(^kerosene 

teT 

For  values  of  dynamic  viscosity  of  =  2.25  x  10"3  Pa-s  and 
rjw  =  1.0  x  10“3  Pa-s,  we  obtain 

3  .0476 


^kerosene 


(6-4) 
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which  is  in  good  agreement  with  the  value  measured  in  the  standard  per¬ 
meability  test. 

Note  that  at  the  highest  apparent  fluid  velocity  attained  in  the  soil 
tests  (32  in/s),  the  equivalent  permeability  k‘y^  from  Equation  3-6  is  about 
15%  of  Darcy's  laminar  permeability,  k.  This  sharp  reduction  in  permeability 
illustrates  the  importance  of  the  influence  of  turbulent  flow  on  fluid  friction. 

Several  attempts  were  made  to  perform  dynamic  tests  on  Enewetak  sand 
samples,  similar  to  those  reported  in  Section  6.2  on  the  flat  duct.  However, 
since  the  apparent  flow  area  c1  the  soil  sample  was  more  than  13  times  that  of 
the  flat  duct,  we  were  not  able  to  generate  sufficient  accelerations  in  the 
pore  fluid  to  make  meaningful  measurements  of  the  influence  of  inertia. 

In  order  to  make  meaningful  inertia  measurements  on  the  soil  samples,  two 
modifications  to  the  test  apparatus  are  required.  First,  a  faster  acting 
valve  would  increase  the  rate  of  application  of  the  pore  pressure  gradient  and 
second,  a  direct  measurement  of  flow  on  the  low  pressure  side  of  the  soil 
sample  could  increase  accuracy  of  the  flow  measurements  by  more  than  an  order 
of  magnitude. 
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Figure  6.1.  Schematic  section  view  of  dynamic  permeability  device. 
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Figure  6.2.  Steady  state  data  replotted  in  order  to  measure  the  coefficient  of  permeability 
and  Ward’s  constant,  flat  duct  flow  test. 
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Figure  6.3a.  Measured  pressure  gradient  and  flow  displacement  time 
histories,  flat  duct  test  J12A8. 


130 


Wi.  Acceleration  fa's)  Normalized  Press.  Grad. 
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Figure  6.3b.  Normalized  fluid  pressure  gradient  and  fluid  motion  time 
histories,  flat  duct  test  J12A8. 
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Figure  6.4a.  Measured  pressure  gradient  and  flow  displacement  time 
histories,  flat  duct  test  F16A8. 
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Figure  6.4b.  Normalized  fluid  pressure  gradient  and  fluid  motion 
time  histories,  flat  duct  test  F16A8. 
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Figure  6.5a.  Measured  pressure  gradient  and  flow  displacement  time 
histories,  flat  duct  test  F1688. 
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Figure  6.5b.  Normalized  fluid  pressure  gradient  and  fluid  motion 
time  histories,  flat  duct  test  P16BG. 
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Figure  6.oa.  Measured  pressure  gradient  and  flow  displacement  time 
histories,  flat  duct  test  F19A8. 
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Figure  6.7a.  Measured  pressure  gradient  and  flow  displacement  time 
histories,  flat  duct  test  F1988. 
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SECTION  7 


TWO-PHASE  RESPONSE  OF  SATURATED  LIMESTONE 


7.1.  INTRODUCTION 

Applied  Research  Associates  is  conducting  a  series  of  laboratory  tests  on 
Indiana  limestone,  under  sponsorship  of  Defense  Nuclear  Agency.  The  testing 
program  was  planned  to  provide  data  to  formulate  and  validate  two-phase  models 
of  the  limestone  response.  In  this  section,  we  review  these  test  data  in  con¬ 
junction  with  the  two-phase  modeling  formulations  we  have  developed  under 
AFOSR  sponsorship.  The  Indiana  limestone  tests  provide  a  data  set  against 
which  to  evaluate  the  two-phase  model  in  a  much  less  porous  material  than  we 
have  previously  had  available. 

In  addition,  a  substantial  portion  of  the  Indiana  limestone  testing  has 
been  performed  along  nonconvent ional  specified  strain  paths.  These  strain 
paths  are,  however,  typical  of  strain  paths  induced  by  explosive  loadings,  and 
the  resulting  stress  paths  provide  test  data  which  will  challenge  modeling 
efforts. 

7.2.  LIMESTONE  PROPERTIES  AND  TESTING  OUTLINE 

The  limestone  selected  for  these  tests  was  quarried  from  the  Salem 
limestone  formation  by  the  Elliott  Stone  Company  of  Bedford,  Indiana.  The 
Salem  limestone  is  a  widely  used  building  stone.  Table  7.1  lists  typical 
index  properties  of  the  Salem  limestone.  The  porosity  of  the  limestone 
averaged  about  13.5%,  which  is  significantly  lower  than  the  36%  porosity  of 
the  typical  Enewetak  limestones  reported  in  our  previous  studies. 

Four  basic  types  of  tests  were  run  on  nominal  2  inch  diameter  by  4  inch 
length  cylindrical  samples  of  the  Salem  limestone.  These  included: 

1.  Hydrostatic 

2.  Uniaxial  strain 

3.  Triaxial  compression  at  constant  confining  pressure 

4.  Specified  strain  paths  following  compression  under  uniaxial  strain. 
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Both  drained  and  undrained  tests  with  pore  pressure  measurements  were  con¬ 
ducted  for  each  test  type. 

The  following  subsections  highlight  the  test  results  and  show  comparisons 
to  two-phase  test  simulations  using  our  two-phase  models. 

7.3.  HYDROSTATIC  TESTS 

A  typical  drained  hydrostatic  compression  test  is  shown  in  Figure  7.1. 
This  is  actually  a  composite  of  several  hydrostatic  tests  to  show  unload- 
reload  response.  The  hydrostat  typifies  the  compression  response  of  the  Salem 
limestone.  There  is  an  initial,  nearly  linear,  portion  of  the  loading  curve 
which  is  essentially  elastic  and  has  a  high  modulus.  Bulk  moduli  of  235  kb 
are  typically  measured  on  the  initial  part  of  the  loading.  At  a  pressure  of 
over  1  kb,  the  cementation  begins  to  break  down  and  the  pore  space  begins  to 
be  crushed  out  during  a  much  softer  portion  of  the  loading  curve.  The  incre¬ 
mental  tangential  bulk  modulus  during  this  portion  of  the  loading  is  nearly  an 
order  of  magnitude  softer  than  the  initial  modulus.  At  pressures  approaching 
4  kb,  the  limestone  begins  to  show  stiffening  as  approximately  half  the  pore 
space  has  been  crushed  out.  Over  the  0  to  4  kb  pressure  range,  the  initial 
unloading  modulus  is  approximately  equal  to  the  initial  loading  modulus,  and 
there  is  a  pronounced  hook,  or  heel,  as  the  pressure  drops  toward  zero. 

A  typical  cyclic  undrained  hydrostatic  loading  of  Salem  limestone  is 
shown  in  Figures  7.2  through  7.4.  The  test  data  are  compared  to  a  numerical 
simulation  computed  using  an  improved  version  of  the  code  NKOCP  described  by 
Kim,  Blouin  and  Timian,  1987.  During  the  past  year,  we  have  modified  NKOCP  to 
include  a  nonlinear  unload  capability  in  addition  to  the  original  bilinear 
unload.  The  NKOCP  simulation  used  the  porous  skeleton  properties  obtained 
from  a  previous  drained  test  as  input.  These  are  combined  with  our  fully 
coupled  two-phase  model  as  described  in  Section  2. 

Figure  7.2  shows  the  total  pressure  as  a  function  of  volume  strain.  The 
initial  undrained  modulus  is  significantly  higher  than  the  drained  modulus, 
having  a  value  of  about  340  kb.  Once  the  sample  is  loaded  beyond  the  pressure 
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at  which  cementation  breaks  down,  there  is  significant  hysteresis  due  to  the 
nonlinear  hysteretic  nature  of  the  skeleton.  The  numerical  NKOCP  simulation 
is  an  excellent  match  to  the  test  data. 

The  effective  stress  response  is  shown  in  Figure  7.3.  Note  that  only 
about  half  the  total  pressure  is  carried  by  the  porous  skeleton.  The  hystere¬ 
tic  nature  of  the  skeleton  is  apparent,  but  liquefaction  does  not  occur  during 
the  unloading  as  it  did  in  the  more  porous  Enewetak  limestone.  However,  we 
would  expect  liquefaction  to  occur  if  the  loading  had  been  carried  to  a 
higher  pressure.  Examination  of  the  drained  response  in  Figure  7.1  indicates 
that  above  art  effective  pressure  of  more  than  1.5  kb  the  unload  becomes  much 
more  hysteretic  which  would  produce  liquefaction  in  an  undrained  test. 

Future  tests  will  be  conducted  at  higher  pressure  to  validate  this  prediction. 
The  agreement  between  the  measured  effective  stress  response  and  that  simu¬ 
lated  by  the  NKOCP  model  is  again  excellent. 

Figure  7.4  shows  the  measured  and  simulated  pore  pressure  response  during 
the  hydrostatic  tests.  Again,  agreement  between  the  two  is  excellent.  There 
is  an  apparent  negative  hysteresis  in  the  pore  pressure  response  which  is  a 
result  of  the  rapid  unloading  of  the  skeleton  and  resulting  rapid  release  of 
compressive  strains  in  the  solid  grains.  This  phenomenon  is  fully  explained 
by  Kim,  Blouin  and  Timian,  1987,  and  highlighted  in  Section  1  of  this  report. 

7.4.  UNIAXIAL  STRAIN  TESTS 

A  typical  drained  uniaxial  strain  test  is  plotted  in  terms  of  axial 
stress  versus  axial  strain  in  Figure  7.5.  The  skeleton  response  is  similar  to 
that  of  the  hydrostatic  drained  test  (Figure  7.1)  but  the  initial  elastic 
constrained  modulus,  M,  is  correspondingly  higher  at  about  420  kb.  The  cemen¬ 
tation  breaks  down  at  an  axial  stress  between  1.5  and  2  kb  and  a  nearly 
linear,  but  softer  portion  of  the  loading  occurs  as  the  porosity  is  crushed 
out  of  the  sample.  The  incremental  tangential  constrained  modulus  along  this 
portion  of  the  loading,  about  50  kb,  is  nearly  an  order  of  magnitude  less  than 
the  initial  modulus.  At  an  axial  stress  of  about  3.5  kb  the  skeleton  begins 
to  stiffen  as  a  significant  volume  of  the  pores  have  been  crushed  out.  The 
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initial  unloading  modulus  tends  to  about  equal  the  initial  elastic  loading 
modulus,  for  peak  axial  stresses  of  less  than  about  4  kb.  At  higher  peak 
stresses,  the  initial  unloading  modulus  appears  to  be  increasing  with 
increasing  peak  axial  stress.  Also  note  that  the  uniaxial  unloading  cannot  be 
maintained  back  to  zero  axial  stress  because  of  plastic  strains  occurring 
during  loading  of  the  sample.  At  some  point  in  the  unloading,  the  axial 
stress  drops  to  the  confining  pressure.  In  order  to  continue  the  unload  and 
maintain  no  lateral  strain,  it  would  be  necessary  to  apply  axial  expansion  to 
the  samples,  a  condition  which  could  not  be  achieved  during  these  tests. 

Thus,  as  noted  in  Figure  7.5,  the  last  portion  of  the  final  unload  is  actually 
hydrostatic,  rather  than  uniaxial. 

The  relationship  between  the  axial  stress  and  the  confining  stress  during 
the  drained  uniaxial  strain  test  is  shown  in  Figure  7.6.  During  the  initial 
elastic  loading  the  apparent  Poisson's  ratio,  ua,  defined  as 

(’-1) 


(T-2) 


where  AK0  a  the  incremental  coefficient  of  lateral  earth  pressure 
&ar  b  incremental  change  in  confining  stress 
Aoa  a  incremental  change  in  axial  stress 

is  equal  to  about  0.25.  During  the  portion  of  the  loading  bey&vtf  the  break¬ 
down  of  cementation,  the  material  responds  more  nearly  hydrostatically,  and 
the  apparent  Poisson's  ratio  increases  to  0.42.  During  the  uniaxial  unload, 
the  apparent  Poisson's  ratio  drops  to  0.32.  The  hydrostatic  portion  of  the 
final  unload  begins  at  a  confining  stress  of  about  2  kb. 


Stress  difference  versus  axial  (or  volume)  strain  is  plotted  for  the 
uniaxial  test  in  Figure  7,7.  At  a  stress  difference  of  about  1  kb,  the 
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cementation  collapses  and  there  is  only  d  slow  further  increase  in  stress 
difference.  In  terms  of  stress  difference,  the  unloadings  reach  a  value  of 
zero  at  which  point  the  sample  is  under  a  hydrostatic  compressive  load,  as  can 
be  seen  in  Figure  7.6. 

The  mean  stress  versus  volume  strain  from  the  uniaxial  test  is  shown  in 
Figure  7.8.  It  is  compared  to  the  drained  hydrostat  from  Figure  7.1.  During 
the  elastic  portion  of  the  loading,  and  at  mean  stress  above  about  3.5  kb,  the 
response  is  essentially  identical  to  the  hydrostat.  However,  the  collapse  of 
the  cementation  occurred  at  consistently  lower  mean  stresses  in  the  uniaxial 
tests  because  the  material  is  subjected  to  an  additional  component  of  shear 
stress.  This  results  in  a  deviation  between  the  two  data  sets  over  the  inter¬ 
mediate  mean  stress  range. 

An  undrained  uniaxial  strain  loading  of  Salem  limestone  is  shown  in 
Figure  7.9  through  7.11.  In  Figure  7.9,  parts  a  through  c,  the  total  axial 
stress,  effective  axial  stress  and  pore,  pressure  are  plotted  against  axial 
strain,  respectively.  A  numerical  simulation  was  generated  using  NKOCP  in 
which  a  set  of  drained  uniaxial  test  data  were  used  as  inputs  to  NKOCP.  The 
agreement  between  the  total  stress  response  of  the  test  and  the  numerical 
simulation  is  excellent.  There  are,  however,  differences  in  the  effective 
stress  response  and  the  pore  pressure  response  between  the  test  and 
simulation.  These  are  primarily  due  to  differences  in  the  skeleton  properties 
of  the  two  samples  being  compared,  as  is  apparent  in  Figure  7.9b.  The  softer 
skeleton  in  the  undrained  test  resulted  in  a  higher  portion  of  the  total 
stress  being  carried  by  the  pore  pressure  in  that  test. 

The  relationship  between  the  axial  and  radial  stress  during  the  undrained 
uniaxial  strain  loading  is  shown  in  Figure  7.10  in  terras  of  both  effective  and 
total  stress.  Because  of  the  stress  component  carried  by  the  pore  fluid,  the 
apparent  values  of  K0,  in  terms  of  total  stress,  are  significantly  higher  than 
values  in  terms  of  effective  stress. 

The  stress  paths  for  the  undrained  uniaxial  loading  are  shown  in  Figure 
7.11  in  terras  of  both  effective  stress  and  total  stress.  The  effective  stress 
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closely  matches  those  from  drained  uniaxial  strain  tests.  It  falls  just  below 
the  drained  ultimate  strength  envelope  obtained  from  triaxial  compression 
tests.  The  total  stress  path  is  much  flatter  than  the  effective  stress  path 
because  of  the  significant  portion  of  total  stress  carried  by  the  pore  fluid. 

7.5.  TRIAXIAL  COMPARISON  TESTS 

Data  from  the  drained  triaxial  compression  tests  on  Salem  limestone  at 
constant  confining  pressure  are  presented  in  Figures  7.12  through  7.15. 

Figure  7.12  shows  stress  difference  versus  axial  and  radial  strain  for  a 
series  of  tests  at  six  different  confining  pressures  ranging  from  0  to  4  kb. 

At  confining  pressures  of  H  kb  and  less  the  samples  fail  in  a  brittle  manner. 
At  higher  confining  pressures  the  samples  fail  in  a  ductile  mode  with  stress 
difference  continuing  to  increase  with  increasing  axial  strain.  The  initial 
modulus  indicated  by  the  slopes  of  the  stress-strain  curves  is  comparatively 
high  for  cafsples  in  the  brittle  failure  regime  at  low  confining  pressures. 
Above  the  point  where  cementation  has  been  crushed  by  the  initial  confining 
pressure,  the  initial  moduli  are  softer  than  those  in  the  brittle  regime.  The 
tests  at  2  and  3  kb  confining  pressure  exhibit  lower  initial  moduli  than  the 
tests  at  lower  confining  pressures.  As  the  initial  confining  pressure  is 
increased  further,  the  modulus  increases  due  to  crush  out  of  the  pore  space  by 
the  hydrostatic  pressure  applied  prior  to  starting  the  triaxial  shearing. 

The  volume  strain  response  during  the  triaxial  shear  portion  of  the 
loading  is  shown  in  Figure  7.13  for  the  six  tests.  In  all  cases,  the  samples 
undergo  compressive  volume  strains  during  the  initial  portion  of  the  shearing. 
As  the  ultimate  failure  stress  is  approached,  the  rock  begins  to  dilate,  as 
indicated  by  the  incremental  volume  expansion.  The  sample  at  2  kb  confining 
pressure  shows  the  greatest  amount  of  compression  during  shear.  At  higher 
confining  pressures  the  compressive  volume  strain  tends  to  decrease,  with  the 
sample  at  4  kb  experiencing  nearly  a  constant  volume  shear. 

The  ultimate  strength  of  the  Salem  limestone  from  the  triaxial  tests  is 
plotted  in  Figure  7.14.  For  purposes  of  this  figure,  ultimate  strength  is 
defined  either  as  the  maximum  stress  difference  reached  during  the  tests  for 
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samples  which  failed  in  the  brittle  regime,  or  as  the  stress  difference  at  15% 
axial  strain  for  samples  at  higher  confining  pressures  in  the  ductile  regime. 
We  have  also  included  test  data  from  Salem  limestone  obtained  by  the  Waterways 
Experiment  Station  (WES)  from  Zelasko,  1988.  The  test  data  define  an  ultimate 
failure  envelope  which  is  linear  below  about  2  kb  stress  difference  and  has  a 
cohesion  of  c  =  0.16  kb  and  a  friction  angle  of  49°  (<f>  =  29°,  in  the  tradi¬ 
tional  Mohr-Coulomb  stress  space).  Above  a  stress  difference  of  2  kb,  the 
envelope  is  concave  downward  toward  the  mean  stress  axis. 

By  plotting  the  stress  difference  at  which  the  incremental  volume  strain 
changes  from  compressive  to  di latent  (i.e.  where  the  slopes  become  infinite  in 
Figure  7.13),  we  can  plot  a  contour  which  defines  the  onset  of  di latency. 

This  is  shown  in  Figure  7.15,  and  lies  just  beneath  the  ultimate  failure  sur¬ 
face  plotted  in  Figure  7.14.  In  other  words,  to  reach  the  ultimate  failure 
surface,  the  rock  must  incrementally  dilate  and  the  incremental  dilation 
begins  when  the  rock  reaches  the  stress  contour  shown  in  Figure  7.15. 

The  undrained  triaxial  compression  test  data  are  presented  in  Figures 
7.16  through  7.19,  The  stress  difference  as  a  function  of  axial  and  radial 
strains  is  shown  for  tests  at  four  confining  pressures  in  Figure  7.15.  All 
samples  exhibit  very  similar  behavior  with  a  rapid  increase  in  stress  dif¬ 
ference  until  failure  at  an  axial  strain  ranging  from  0.5  to  1.0%,  The 
strength  is  only  slightly  dependent  on  the  total  confining  pressure.  There  is 
no  further  increase  in  stress  difference  with  increasing  axial  strain. 

The  effective  stress  paths  for  the  four  undrained  tests  are  shown  in 
Figure  7.17.  The  initial  effective  stresses  following  the  hydrostatic  portion 
of  the  test  are  approximately  half  of  the  total  confining  pressure.  The  maxi¬ 
mum  effective  stress  difference  in  each  test  agrees  reasonably  well  with  the 
drained  ultimate  strength  envelope  from  Figure  7.14.  In  terms  of  total 
stress,  the  strength  envelope  would  exhibit  a  friction  angle  of  only  about 
10®,  as  indicated  in  Figure  7,17  (vs  49®  for  the  drained  envelope).  Note  that 
the  effective  stress  paths  tend  to  be  concave  toward  the  mean  effective  stress 
axis  at  low  confining  pressures,  but  concave  toward  the  stress  difference  axis 
at  higher  confining  pressures.  This  difference  in  stress  indicates  that  pore 
pressures  are  increasing  more  rapidly  in  the  higher  confining  pressure  tests. 
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Stress  difference  as  a  function  of  volume  strain  for  the  undrained  tests 
is  shown  in  Figure  7.18.  As  would  be  expected,  there  is  a  slightly 
compressive  volume  strain  occurring  during  the  initial  portion  of  the  loading. 
Once  the  effective  stresses  reach  the  strength  envelope,  however,  there  is  a 
sudden  failure  and  the  samples  begin  to  dilate  strongly.  Dilation  occurs 
throughout  the  remainder  of  the  test,  with  no  further  increase  in  stress  dif¬ 
ference  . 

A  corresponding  plot  of  pore  pressure  as  a  function  of  volume  strain  is 
shown  in  Figure  7.19.  There  is  an  increase  in  pore  pressure  during  the  ini¬ 
tial  loading,  followed  by  a  sudden  leveling  off  once  failure  occurs.  During 
the  di latent  portion  of  the  test  the  pore  pressure  remains  nearly  constant  for 
the  tests  at  the  higher  confining  pressures  and  decreases  slowly  for  the 
0.5  kb  confining  pressure  test. 

The  behavior  demonstrated  by  the  test  data  in  Figures  7.18  and  7.19  and 
contradicts  what  we  would  expect  from  effective  stress  theory  and  demonstrates 
a  shortcoming  of  undrained  triaxial  strength  tests.  The  volume  strain  is 
measured  at  the  midpoint  of  the  cylindrical  test  specimen.  At  the  midpoint, 
where  the  shear  failure  is  occurring,  the  sample  is  dilating  rapidly.  Local 
pore  pressures  within  this  shear  zone  should  be  dropping  rapidly  as  volume 
expansion  occurs.  As  pore  pressures  drop,  effective  stresses  and  the  stress 
difference  should  both  increase.  Figures  7.16  through  7.19  clearly 
demonstrate  that  this  does  not  occur  during  the  latter  stages  of  the  undrained 
tests.  We  believe  the  reason  that  pore  pressure  doesn't  decrease  is  that 
total  volume  strain  of  the  specimen  remains  nearly  constant  during  shearing. 
For  the  total  volume  strain  to  remain  constant,  the  ends  of  the  sample, 
constrained  radically  by  the  steel  end  caps,  must  be  undergoing  volumetric 
compaction  while  the  center  is  experiencing  the  measured  volumetric  expansion. 

In  essence,  we  believe  that  the  radial  constraint  imposed  by  the  sample 
end  caps  prevents  the  sample  strength  from  increasing  during  the  later  stages 
of  shear.  This  is  a  seriou*.  shortcoming  of  the  undrained  triaxial  test  which 
underestimates  the  sample  strength  at  higher  shear  strains  and  could  lead  to 
gross  errors  in  predicting  the  late  time  in  situ  response  of  saturated  porous 
media. 
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7.6.  STRAIN  PATH  TESTS 


A  series  of  drained  and  undrained  load-unload  cycles  was  performed  on 
Salem  limestone  following  predesignated  strain  paths.  The  strain  paths  were 
chosen  to  be  representative  of  strain  paths  from  various  depths  beneath 
shallow  buried  explosive  detonations  as  predicted  by  finite  difference  calcu¬ 
lations  performed  by  several  of  the  national  laboratories.  Figure  7.20  is  a 
schematic  view  of  the  strain  paths.  All  paths  undergo  an  initial  compression 
under  uniaxial  strain  conditions,  then  from  a  designated  point  on  the  uniaxial 
loading  curve,  one  of  four  paths  is  followed: 

1.  Path  A:  a  constant  volume  path  in  which  axial  strain  continues 

in  compression  while  radial  strain  is  allowed  to  expand  at  half  the 
rate  of  axial  compression  resulting  in  a  constant  volume  shearing; 

2.  Path  8:  a  constant  axial  strain  path  in  which  axial  strain  is  held 
constant  while  radial  expands; 

3.  Path  C:  an  iso-unloading  path  in  which  both  axial  and  radial  strains 
expand  at  the  same  rate; NS 

A.  Path  0:  a  uniaxial  strain  unload  which  is  the  same  as  that  used  in 
the  uniaxial  strain  tests  described  in  Section  7.4. 

The  actual  drained  strain  paths  and  resultant  stress  Paths  for  paths  A, 

8,  and  G  are  plotted  in  Figures  7.2la,  7.22a  and  7.23a,  respectively.  The 
Type  A  constant  volume  stress  paths,  sMwn  in  Figure  7.21b,  show  increasing 
shear  stress  with  decreasing  mean  stress,  approaching  the  constant  volume 
strength  contour  (from  Figure  7.15)  along  a  path  nearly  perpendicular  to  it. 
Similar  tests  run  by  WES  to  larger  shear  strains  showed  some  further  gain  in 
shear  strength  as  the  stress  path  turned  upward  and  toward  the  right  as  it 
moved  along  the  constant  volume  contour. 

The  Type  8,  constant  axial  strain  stress  paths,  are  shown  in  Figure 
7.22b.  These  tend  to  move  laterally  toward  the  ultimate  strength  envelope 
without  further  increase  in  stress  difference.  In  all  cases,  the  stress  paths 
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extend  slightly  beyond  the  strength  envelope,  which  indicates  that  the 
material  has  been  strengthened  somewhat  by  the  compaction  process  during 
uniaxial  strain  loading. 

The  Type  C,  iso-unloading  stress  paths,  are  shown  in  Figure  7.23b.  These 
tend  to  unload  down  the  uniaxial  loading  stress  paths  without  the  sample 
reaching  failure. 

Figure  7.24  is  an  idealized  schematic  representation  of  the  stress  paths 
resulting  from  the  four  prescribed  strain  paths.  The  stress  paths  represent  a 
range  of  material  response  which  will  be  difficult  to  match  with  most  existing 
models. 

Undrained  strain  paths  for  samples  loaded  along  Path  B  (constant  axial 
strain)  and  Path  C  (iso-unload)  are  shown  in  Figures  7.25  and  7.26.  In  the 
Stress  path  for  the  Type  B  loading  of  Figure  7.25b,  the  effective  stress 
response  is  very  similar  to  that  of  the  drained  tests  in  Figure  7.22b.  Note 
that  both  the  effective  stress  path  and  the  total  stress  path  converge  on  the 
strength  envelope  indicating  that  pore  pressure  has  dropped  to  zero  at  this 
point.  Thus,  the  drained  strength  envelope  serves  both  the  effective  stress 
and  total  stress  response  for  the  Type  B  strain  path. 

The  effective  stress  response  for  the  Type  C  strain  path,  shown  in 
Figure  7.26b,  is  also  very  similar  to  that  of  the  drained  tests  shown  in 
Figure  7.23b.  Note  that  pore  pressure  also  drops  to  zero  for  this  load  path 
prior  to  complete  unloading  and  that  the  sample  does  not  undergo  failure. 

Figure  7.27  shows  hypothesized  strength  envelopes  in  terms  of  both  effec¬ 
tive  and  total  stresses  for  undrained  Type  A  and  B  strain  paths  and  triaxial 
compression  tests.  These  were  constructed  using  NKOCP  numerical  simulations 
of  the  initial  hydrostatic  and  uniaxial  strain  portions  of  the  loadings,  then 
approximating  the  shear  loadings  using  effective  stress  theory  and  adaptation 
of  our  two-phase  model.  This  figure  demonstrates  that: 

1.  strength  envelopes  expressed  in  terms  of  effective  stress  do  not  vary 
with  the  strain  path  and  match  the  envelopes  determined  from  drained 
tests;  and 
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2.  strength  envelopes  expressed  in  terms  of  total  stress  vary  dramati¬ 
cally,  being  dependent  on  the  strain  path. 

Because  of  the  dependence  of  undrained  strength  on  the  strain  paths,  the 
appropriateness  of  using  equivalent  single-phase  models  for  ground  shock 
calculations  in  saturated  media  (a  common  practice)  is  questioned.  Such 
models  cannot  make  direct  compensation  for  the  influence  of  various  strain 
paths  on  strength  and  post  failure  response. 
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Table  7.1.  Typical  Salem  Limestone  Properties. 


Salem  Limestone 
(Avq.  6  samples) 


Ory  Density  (lb/ft3)  (kg/m3) 

147.4 

2359 

Porosity 

.135 

Specific  Gravity  of  Grains 

2.728 

Permeability  (cm/s) 

5.3  X  10-5 
(1  sample) 

Ultrasonic  P  Wave  Velocity  (ft/s)  (m/s) 

14806 

4514 

Yield  Stress  (unconfined)  (lb/in2)  (kb)-. 

5540 

.382 

Ultimate  Stress  (unconfined)  (lb/in2)  (kb) 

7752 

.535 

Young's  Modulus  (lb/in2)  (kb) 

3.9  X  106 

270 

Yield  Strain  (unconfined)  (Axial,  %) 

,168 

Ultimate  Strain  (unconfined)  (Axial,  %) 

.261 
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al  drained  hydrostat  on  Salem  limestone 


Figure  7.2.  Comparison  of  total  stress  from  an  undrained  hydrostatic  compression  test  with 
a  two-phase  simulation  using  skeleton  properties  from  drained  test  data. 
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Figure  7.3.  Comparison  of  effective  stress  from  an  undrained  hydrostatic  compression 
test  with  a  two-phase  simulation  using  skeleton  properties  from  drained 
test  data. 
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Figure  7.4.  Comparison  of  pore  pressure  from  an  undrained  hydrostatic  compression  v 
with  a  two-phase  simulation  using  skeleton  properties  from  drained  test 
data. 


2350  kg/m 


Figure  7.5.  Typical  drained  uniaxial  strain  loading  on  Salem  limestone. 
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Figure  7.6.  Drained  uniaxial  strain  test  on  Salem  limestone;  relationship  between 


difference  and  volume  strain. 


test  on  Salem  limestone;  relationship  between  true  mean  stress 


o 
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Figure  7.3b.  Comparison  of  effective  stress  from  an  undrained  uniaxial  strain  test  with 
a  two-phase  simulation  using  skeleton  properties  of  Salem  limestone  from 
drained  test  data. 


Figure  7.9c.  Comparison  of  pore  pressure  from  an  undrained  uniaxial  strain  test  with 
a  two-phase  simulation  using  skeleton  properties  of  Salem  limestone  from 
drained  test  data. 


Radial  Stress  (Kbar) 


Figure  7.11.  Stress  paths  for  undrained  uniaxial  strain  loading  of  Salem  limestone  -  Test  L11A8 


Figure  7.12.  Ax  id  1 -radia ?  strain  response  cf  drained  Salem  limestone  in  triaxial  compression 
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Figure  7.13.  Relationship  between  stress  difference  and  volume  strain;  Salem 
limestone  in  drained  triaxial  compression. 
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Figure  7.15.  Contour  for  the  onset  of  di latency  in  drained  triaxial  compression  of  Salem 
1 imestone. 
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Figure  7.16.  Axial  and  radial  Strain  respo; 
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Figure  7.17.  Effective  stress  paths  for  Salem  limestone  in  undrained  triaxial  compression 
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Figure  7.20.  Strain  paths  specified  for  Salem  limestone  tests. 


Figure  7.21a.  Type  A  strain  paths  for  drained  Salem  limestone. 


Strength  Envelope 


Figure  7.21b.  Stress  paths  for  drained  Type  A  strain  paths  on  Salem  limestone. 
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Figure  7.22b.  Stress  paths  for  drained  Type  B  strain  paths  on  Salem  limestone. 
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Figure  7.23a.  Type  C  strain  paths  for  drained  Salem  limestone. 
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Figure  7.23b.  Stress  paths  for  drained  Type  C  strain  paths  on  Salem  limestone. 
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Figure  7.25b.  Stress  path  for  undrained  Type  B  strain  path  test  for  Salem  limestone. 
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Figure  7.26b.  Stress  patfts  for  undrained  Type  C  strain  path  test  for  Salem  limestone. 


Drained  Ultimate  Strength  J  Triaxia! 

Undrained  Ultimate  Effective  Strength  (  Shear 


Figure  7.27.  Hypothesised  strength  envelopes  for  Type  ft  and  Type  B  strain  paths  on  Salem  limestone 
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IBATCH  (15) 


iBATCH  =  0  Terminal  interactive  job 
1  Batch  job 

CARD  2 

Main  Title  (up  to  80  characters) 

CARO  3 


Subtitle  (up  to  80  characters) 


CARO  4  General  Options 


NF,  NTCSF,  ISFG,  IP,  NLNR,  NFG,  ICOPL  (715) 

Options  in  Dynamic  Analysis 

Nf-  3  1  One-phase  solid  dynamic  analysis 

2  One-phase  fluid  dynamic  analysis 

3  Multiphase  dynamic  analysis 

Analysis  Type 

NTCSF  =  1  Static  analysis 

2  Consolidation  analysis 

3  Dynamic  analysis 

Initial  Stress  Conditions 
ISF6  »  0  No  initial  stress 

1  Specified  effective  stress  and  pressure 

2  Specified  effective  stress 

3  Specified  pore  pressure 

4  Imposed  excess  pore  pressure 

Stress  and  Strain  Conditions 
IP  »  1  Plane  stress 

2  Plane  strain 

3  Axial  symmetry 

4  Spherical  symmetry 

Skeleton  Material  Model 
NLNR  &  0  Linear  elastic  material  model 

1  Decoupled  elesto-plastic  model 

4  Uniaxial  strain  model 

5  Generalized  Hoek  and  3rown  rock  model 

6  ARA  3- Invariant  model 

Loading  Type 

NFG  “  1  Specified  base  accelerations  (not  available) 

2  Specified  pressure  time  history 

3  Specified  velocity  time  history  (not  available) 

4  Both  pressure  and  velocity  time  histories  (not  available) 

Fluid  -  Grain  -  Skeleton  Compressibility  model 
ICOPL  n  0  Decoupled  model 

1  Fully  coupled  model 


CARO  5  Global  Calculator  Parameters 

NCYCL ,  DT,  NUPDAT,  ITER  (15,  F10.0,  215) 

NCYCL :  Number  of  cycles  (total  number  of  time  steps) 

DT:  Global  time  step  (duration  of  each  cycle) 

NUPDAT:  Number  of  cycles  between  updates  to  global  stiffness  matrix 
(ITER  =  0  for  NUPDAT  >  1) 

ITER:  Number  of  updates  of  global  stiffness  matrix  within  each 

cycle  (NUPDAT  =  1) 


CARD  6  Mesh  and  Material  Parameters 

NoMNf- ,  NUMEL,  NUMMAT,  NHSIZE,  NPSIZE  (515) 

NUMNP :  Number  of  nodes 

NUMEL:  Number  of  elements 

NUMMAT:  Number  of  different  materials 

NHSIZE:  Maximum  number  of  histeretic  variables  (for  ARA  3-Invariant 
model,  use  NHSIZE  =  14) 

NPSIZE:  Maximum  number  of  material  parameters  (for  ARA  3-Invariant 
model,  use  NPSIZE  »  100) 


CARO  1  Pressure  Loading  Functions 

NUMLP,  NUMLH,  NUMTP,  NTYPE,  OTXX  (418,  F10.0) 

NUMLP:  Total  number  of  degrees  of  freedom  at  which  input  pressure 

time  history  are  specified 
NUMLH:  Number  of  input  pressure  time  histories 

NUMTP:  Number  of  pressure  time  pairs  in  every  input  pressure  time 

history 

NTYPE  a  0;  Constant  time  increments  in  pressure  time  history 
it  Specified  times  in  pressure  time  history 
DTXX:  Constant  time  increment  in  the  input  pressure  time 

history  (NTYPE  ■  0) 
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CARO  8 


Output  Stress  and/or  Motion  Profile  Specifications 

A.  NPFL,  NDC,  NSG,  NPRINT,  NPEL,  NPMT  (615) 


NPFL  =  0  Both  motions  at  all  nodes  and  stresses  at  all  elements 

1  Motions  and  stresses  at  specified  nodes  and  elements 
respectively 

2  Motions  at  all  nodes 

3  Stresses  at  specified  elements 

4  Stresses  at  all  elements 

NDC  =  0  Write  stress/displacement  profile  output  to  hard  disc 

1  Write  stress/displacement  profile  output  to  floppy  disc 


NSG 


0  Write  stress/displacement  profile  output  in  one  file 

1  Write  stress/displacement  profile  output  in  specified 

files 


NPRINT:  Number  of  cycles  between  each  output  profile 

NPEL:  Number  of  elements  in  output  stress  profile 

NPMT:  Number  of  nodes  in  output  motion  profile 

NOTE:  If  NDC  a  l  and  NSG  =  0,  the  program  writes  displacement/stress 
profiles  on  hard  disk  under  the  name  "DISTR" 

B.  If  NPFL  a  i  or  3,  otherwise  go  to  Card  8C. 

NPRT(I),  1=1,  NPEL  (free  format) 

NPRT(I);  List  specified  element  numbers  in  sequential  order 

C.  If  NPFL  =  1,  otherwise  go  to  next  card 
NMP(I),  1=1,  NPMT  (free  format) 

NPH(I):  List  specified  node  numbers  in  sequential  order 
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CARO  9 


Output  Stress  and/or  Motion  Time  History  Specifications 


A.  NTHS,  NHPEL,  NHPMT  (315) 

Time  history  options 

NTHS  =0  Do  not  print  time  history  data,  go  to  next  card 

1  Motion  time  histories 

2  Stress  time  histories 

3  Both  motion  and  stress  time  histories 

NHPEL:  Number  of  elements  at  which  stress  time  histories  are 

required 

NHPMT:  Number  of  nodes  at  which  motion  time  histories  are 
required 

B.  If  NTHS  a  2  or  3,  otherwise  go  to  9C. 

NHPRT(I),  I=l,  NHPEL  (free  format) 

NHPRT(I):  List  specified  element  numbers  in  sequential  order 

C.  If  NTHS  »  1  or  3,  otherwise,  go  to  next  card 
NHPM(I),  1=1,  NHPMT  (free  format) 

NHPM( I ) :  List  specified  node  numbers  in  sequential  order 
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CARO  10  Numerical  Time- Integration  Options 


TETA,  BETA,  GAMA,  ALPA,  (4F10.0) 
TETA:  0 

BETA:  0  Refer  to  Table  1. 
GAMA:  y 
ALPA:  a 


TABLE  1.  Values  of  0  and  0  for  y  = 

Integration  Method 
(1) 

1/2*  (a  = 

0 

(2) 

0  by  Default) 

e 

(3) 

Explicit  second  central  difference 

0 

1.0 

Fox -Goodwin 

1/12 

1.0 

Linear  Acceleration 

1/6 

1.0 

Newmark's  constant  acceleration 

1/4 

1.0 

Wilson 

1/6 

2.0 

Stiff  linear  acceleration 

1/6 

1.5 

*y  =>  1/2  indicates  no  damping 

y  >  1/2  introduces  numerical  damping  and  <5  *  (y  ♦  l/2)2/4 
a  Method  available  soon 

For  more  information,  see:  Ghaboussi  and  Wilson,  "Variational 
Formulation  of  Oyanmic  of  Fluid  Saturated  Porous  Elastic  Solids," 
ASCE  Engineering  Mechanics  Journal,  August  1972. 
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CARO  11  Nodal  Coordinates  and  Degree  of  Freedom  Specifications 
For  each  node: 

NOOE,  ISX,  ISY,  IFP,  XA,  YA  (415,  2F10.0) 


Y 


NODE :  Node  Number 

ISX:  Specifies  skeleton  X  Degrees  of  Freedom  (DOF) 

ISY:  Specifies  skeleton  Y  DOF 

IFP:  Specifies  DOF  for  pore  fluid  pressure 

ISX,  ISY  ®  0  Free  to  move  in  specified  direction 
1  Fixed  in  specified  direction 
IFP  a  0  Unknown  pore  fluid  pressure 
1  Zero  pore  fluid  pressure 

XA:  X  Coordinate 

YA:  V  Coordinate:  Note  for  IP  =  4  (l-D  spherical  analysis)  set  the 

mesh  height  equal  to  1.0. 
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CARO  12  Element  Specifications 


For  each  element; 

NEL,  MAT,  KS,  KF,  INTR,  INTS,  I,  J,  K,  L,  Ml,  M2,  M3,  M4  (1415) 


NEL:  Element  number 

HAT:  Material  property  number 

KS  w  0  Element  has  solid  phase 

l  Element  has  no  solid  phase 

KT  *  o  Element  has  fluid  phase 

'V-  Element  has  no  fluid  phase 
INTR:  Number  of  integration  points  in  r-direction 

INTS:  Number  of  integration  points  in  s-direction 

Node  numbers  at  element  corner 
M1,M2,H3.M4:  Node  numbers  at  element  midside 
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CARO  13  MATERIAL  PROPERTIES 


A.  Title  (up  to  80  characters) 

B.  General  Property 
POR,  GF  (Free  Format) 

POR:  Initial  porosity 
GF:  Unit  weight  of  pore  fluid 

C.  Permeability  Property 

NP,  RK,  BETHF  (free  format) 

NP  =  0:  Constant  permeability 
=  1:  Nonlinear  permeability 
RK  :  Oarcy's  coefficient  of  permeability 
BETHF:  Ward's  fluid  friction  coefficient  for  turbulent  flow 


0.  Grain  Property 

NG,  BKG,  ROG,  CO,  VO,  S,  PB  (free  format) 


NGafl 

1 

BKG 

ROG 

CO 

VO 

S 


P0: 


Constant  grain  •, modulus 

Nonlinear  grain  modulus 

Initial  bulk  modulus  of  grain 

Initial  mass  density  of  grain 

Initial  wave  velocity  at  relatively  low  pressure 

Initial  Poisson's  ratio 

Experimentally  determined  constant  relating 
loading  wave  velocity  to  peak  particle  velocity 
(generally  about  equal  to  1.5  for  most  rocks  and 
minerals) 

Threshold  pressure  beyond  which  materials  tend  to 
behave  like  fluids 


E.  Fluid  Property 

NW,  8KW,  ROF,  SO,  HC,  PAO,  T  (free  format) 

NW  a  0:  Constant  fluid  modulus 
I:  Nonlinear  modulus 
BKW:  Initial  bulk  modulus  of  pore  fluid 
ROF:  Initial  mass  density  of  pore  fluid 
SO:  Initial  degree  of  saturation 
HC:  Coefficient  of  solubility  (Henry's  constant) 

PAO:  Initial  pore  air  pressure  (absolute) 

T;  Pressure  difference  between  the  air  and  pore  water  due  to 
surface  tension 
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F.  Skeleton  Property 


a.  General  Property 

NS,  BKS,  GS,  TENS,  STIFAC,  SHEFAC,  PMN  (free  format) 

NS  =  0:  Linear  elastic  modulus 
1:  Nonlinear  modulus 
BKS:  Initial  bulk  modulus  of  skeleton 
GS:  Initial  shear  modulus  of  skeleton 
TENS  =  0.0:  No  tension  cut-off 
1.0:  Tension  cut-off 

STIFAC:  Factor  which  reduces  normal  stiffness  once  tensile 
strength  is  exceeded  (example:  reduced  modulus  » 
original  modulus/STIFAC) 

SHEFAC:  Factor  which  reduces  shear  stiffness  once  shear 
strength  is  exceeded  (see  example  above) 

PMN:  Tensile  strength  (tensile  stress  is  positive) 

b.  Material  Model  Specification 
NQ0N0  (free  format) 


HOONO  a  0: 

1: 

4; 

5: 

6: 


Linear  elastic  material  model 
Decoupled  elasto-plastic  model 
Uniaxial  strain  model 
Generalized  Hoek  and  Brown  rock  model 
AftA  3 -Invariant  mode! 


c.  Decoupled  Elasto-Plastic  Model  (OCOUP) 

If  KQONO  is  not  equal  to  1,  go  to  next  card  group. 
1.  Shear  Stress-Strain  Behavior 


GP,  A,  8  (free  format) 

GP:  Plastic  shear  modulus 
A:  a  (constant  defining  proporitonal  elastic  limit) 

B:  b  (fraction  of  plastic  shear  modulus  afte  failure) 
(refer  to  Model  OCOUP  by  Kim,  Qlouin  and  Tioian, 
1986) 
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fc'-anSufcrt 


**0ESCRI8E  VOLUMETRIC  BEHAVIOR** 


2.  NLPC  (15) 

3*  P^»  > 

p2>  bL2 

“  ”  >  "NLPC"  cards  with  format  (2F10.0) 

I  for  each  card 

pn  &Ln  / 


4.  NUPC  (15) 

5.  Pj,  Bui 

p2*  3U2 

«■* 

pn*  BUn 


“NUPC"  cards  with  format  (2P10.0) 
for  each  card 


* 


a 

i 
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NLPC:  Number  of  pressur^/modulus  pairs  describing 
the  loading  bulk  moauiys 
BLi :  Loading  bulk  modulus  at  mean  pressure 

NUPC :  Number  of  pressure/modulus  par'is  describing 
*he  unloading  bulk  modulus 
Bg^ :  Unloading  bulk  modulus  at  mean  pressure  P, 
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**OESCRIBE  DEVIATORIC  BEHAVIOR** 


6.  NTPC  (15) 

7.  Pi,  Tyi i  TF1 
p2'  TY2'  Tf 2 

pn>  pYn>  pFn 


\ 


|  “NTPC"  cards  with  format  (3F10.0)  for 

>  each  card 


NTPC:  Number  of  presure/yield  stress/failure  stress 
points  describing  the  yield  and  failure 
envelopes  as  a  function  of  mean  pressure 

Ty^  Tpi :  Octahedral  shear  stress  at  yield!  and 
failure,  respectively,  at  mean  pressure 
Pi 


Uniaxial  Strain  Model  (UNIAX)  Reference:  Kim  and  Blouin,  1986 


If  KODNO  is  not  equal  to  4  go  to  next  card  group 

1.  POSNR,  EQNO,  C,  0,  SVMALL  (free  format) 

POSNR :  Poisson's  ratio 

EQNO  =  Is  Unloading  modulus  as  a  function  of  effective 
vertical  stress 

2:  Unloading  modulus  as  a  function  of  previous 
maximum  effective  vertical  stress  (<7v'max) 

3:  Unloading  modulus  computed  by  Mu  =  C  (av'max)° 
C,0:  Material  parameters  defining  unloading  modulus 
when  EQNO  -  3 


2.  NLPC  (15) 


3.  S1#  M^ 


Sn,  Mi_n 

4.  NUPC  (15) 

5.  Sj,  Mui 


sn*  ^Un 


"NLPC"  cards  with  format  (2F10.0)  for 
each  card 


'•NUPC"  cards  with  format  (2F10.0)  for 
each  card 


NLPC:  Number  of  vertical  stress- loading  constrained 
modulus  pairs 


Loading  constrained  modulus  at  vertical  stress 

NUPC:  Number  of  vertical  stress-unloading  constrained 
modulus  pairs 


My^:  Unloading  constrained  modulus  at  vertical 

stress  S* 


e.  Generalized  Hoek  and  Brown  Rock  Model 


If  NQDNO  is  not  equal  to  5,  go  to  next  card  group. 


RN,  ALPHA,  BETA,  RLKAPA,  RLK  (free  format) 
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RN:  n 
ALPHA:  a 
BETA:  0 
RLKAPA:  K 
RLK:  K 

Note  chat  the  strength  parameters,  n,  a,  0,  tc,  and  K  are 
described  in  detail  by  Kim  and  Blouin,  1986. 

f.  ARA  3  -  Invariant  Model,  Reference:  Merkle  and  Oass,  1985) 

If  MODNO  is  not  equal  to  6,  go  to  next  card  group. 

1.  APEX,  ATMO  (free  format) 

APEX:  Tensile  strength  measured  along  the  octahedral 
normal  sterss  axis 

ATMO:  Atmospheric  pressure  (0.1013  MPa) 


2.  AKUR,  AN,  APOI  (free  format)  Elastic  Constants 

AKUR:  Constant  for  Young's  modulus  (Kur 
AN:  Constant  .for  Young's  modulus  (n) 

APOI;  (Poisson's  ratio  (u) 

3.  AR,  ACRV  (free  format)  Compressive  Yielding  and 
Hardening  Parameters 

AR:  Cap  axes  ratio 
ACRV;  Number  of  c  and  p  data  set 

4.  AACC(l),  AAPC(l) 

AACC ( 2 ) ,  AAPC(2)  Free  Format  "ACRV"  cards 

»  * 


I 


AACC:  Hardening  constant  (c) 

AAPC;  Hardening  constant  (p) 

6.  AEY,  AMY,  AETA1  (free  format)  Expansive  Yielding 
and  Hardening  Parameters 


AEY:  Yield  constant  (E) 
AMY:  Yield  exponent  (m) 
AETAi:  Failure  constant  (nj) 


6.  ATO,  ARC,  AS6  (free  format)  Plastic  Potential  Constant 
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ATG:  t 
ARG:  R 
ASG:  S 

1.  APBAR ,  AL  (free  format)  Work  Hardening  Constant 

AP8AR :  p 
AL:  J 

Note  that  the  peak  expansive  plastic  work,  Wp,peak, 
is  related  to  initial  confining  stress  as 

Wp,peak  =  p  •  Pa  (p—)* 


8.  AALPH,  ABETA  (fne  format) 

AALPH:  Work  hardening  constant  (a) 

ABETA:  Work  hardening  constant  (0) 

No*e  -'.hat  the  hardening  parameter,  q,  is  related  to  the 
initial  confining  stress  as 


q  o  o  +  6  (£|£) 

and  the  program  MPDAP  computes  the  hardening  parameters 
a  and  b  as  follows: 


a 


a 


*1 


and 


b  a 


l 


q  •  WD,peak 
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CARO  14  Initial  In  Situ  Effective  Stress  Conditions 


If  ISFG  from  Card  4  equals  1  or  2,  otherwise  go  to  next  card 

For  each  element 

SXX,  SYY,  S 17,  SXY  (4F10.0) 

Initial  effectie  skeleton  stresses 

SXX  =  ax'  (normal  stress  in  x  direction) 

SYY  =  Oy'  (normal  stress  in  y  direction) 

S2Z  =  oz'  (normal  stress  in  z  direction) 

SXY  =  rXy  (shear  stress  in  xy  plane) 

CARO  15  Initial  Pore  Fluid  Pressures 

If  ISFG  from  Card  4  equals  1,  3,  or  4,  otherwise  go  to  next  card 
PRF(I)  I  s  l,  NUHEL  (F10.0) 

PRF:  List  of  initial  pore  fluid  pressures  in  each  element,  specified 
sequentially  from  1  to  NUHEL. 

Note  that  tensile  stresses  are  positive. 
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CARD  16 


Input  Loading:  For  specified  pressure  and  flow  time  histories 
(Card  4,  NFG  =  2) 


A.  Header  (80  characters) 

B.  For  each  of  the  NUMLP  nodes  (Card  7)  at  which  a  load  is  applied: 
NODE,  I DOF,  LHNO,  CINT  (3I5,E10.3) 

C.  TD(I ) ,  I  =  1,  NUMTP  (7E10.0)  for  (NTYPE  =  1) 


D.  For  each  of  the  NUMLH  (Card  7}  loading  time  histories: 
DYL(I),  I  =  1,  NUMTP  (7E10.0) 

E.  Comment  Card  (80  characters)  -  to  separate  loadings 


NODE :  Node  number 

IOQF  =>  1  Total  force  acting  on  a  given  node  in  the  x  direction 

2  Total  force  acting  on  a  given  node  in  the  y  direction 

3  Flow  at  the  given  node 

LHNO:  Load  history  number  (one  of  NUMLH  load  histories) 

CINT:  Load  intensity  factor  (used  to  convert  pressure  or  flow  on  a 

given  node,  based  on  contributing  area) 

TD( I ) :  Set  of  specified  times  used  in  all  the  pressure  loading 

histories 

OYL(I):  A  set  of  pressure/flow  magnitudes  input  for  each  loading 

history  at  corresponding  times  T0(I) 
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APPENDIX  B 

STEADY  STATE  FLOW  DATA  THROUGH 
ENEWETAK  BEACH  SAND 
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/  in  ) 


U21C8 


Gradient  (psi  /  in  ) 


/  in  ) 


U24B8 


0  30  eo  80  ISO  ISO 


Time  (s) 
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Apparent  Flow  Displacement  (in) 


U27D8 


Time  (s) 


Time  ($) 


! 

[ 

i 


Apparent  Flow  Displacement  (in)  Pressure  Gradient  (ps.  /  in  ) 


U27B8 


1 


Time  (s) 
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Apparent  Flow  Displacement  (in)  Pressure  Gradient  (psi  /  in  ) 


G9F8 


Time  (s) 
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U27J8 


Time  (s) 


Time  (s) 
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/  in  ) 


U27K8 


/  in  ) 


U27L8 


Gradient  (psi  /  in  ) 


U20B8 


Gradient  (psi  /  in  ) 


U20C8 


Gradient  (psi  /  in  ) 


U20D8 


Apparent  Flow  Displacement  (in)  Pressure  Gradient  (psi  /  in  ) 


U21A8 


Time  (s) 


Time  (s) 
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Gradient  (psi  /  in  ) 


U21B8 


psi  /  in  ) 


U30I8 


U23D8 


psi  /  in  ) 


U27H8 


U30A8 


Apparent  Flow  Displacement  (in)  Pressure  Gradient  (psi  /  in  ) 


U30E8 


Time  (s) 


Time  (s) 
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psi  /  in  ) 


Apparen 


U30H8 


Time  (s) 


242 


